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Abstract 


An  Exact  Ceiling  Point  Algorithm 
for  General  Integer  Linear  Programming 

\  Robert  M.  Saltzman  and  Frederick  S.  Hillier 

Stanford  University,  1988 

This  report  describes  an  exact  algorithm  for  the  pure,  general  integer  linear  program¬ 
ming  problem  ( ILP ).  Common  applications  of  this  model  occur  in  capital  budgeting 
(project  selection),  resource  allocation  and  fixed-charge  (plant  location)  problems.  The 
central  theme  of  our  algorithm  is  to  enumerate  a  subset  of  all  solutions  called  ^feasible 
1-ceiling  points.  A  feasible  1-ceiling  point  may  be  thought  of  as  an  integer  solution  lying 
on  or  near  the  boundary  of  the  feasible  region  for  the  LP-relaxation  associated  with  (ILP). 
Precise  definitions  of  1-ceiling  points  and  the  role  they  play  in  an  integer  linear  program 
are  presented  in  a  recent  report  by  the  authors.  One  key  theorem  therein  demonstrates 
that  all  optimal  solutions  for  an  (ILP)  whose  feasible  region  is  non-empty  and  bounded 
are  feasible  1-ceiling  points.  Consequently,  such  a  problem  may  be  solved  by  enumerating 
just  its  feasible  1-ceiling  points.  Our  approach  is  to  implicitly  enumerate  1-ceiling  points 
with  respect  to  one  constraint  at  a  time  while  simultaneously  considering  feasibility.  Com¬ 
putational  results  from  applying  this  incumbent-improving  Exact  Ceiling  Point  Algorithm 
to  48  test  problems  taken  from  the  literature  indicate  that  this  enumeration  scheme  may 
hold  potential  as  a  practical  approach  for  solving  problems  with  certain  types  of  structure. 


Subject  Classification  for  OR/MS  Index:  Programming  -  Integer  Algorithms; 
Branch-and-bound 


Key  Words:  integer  linear  programming;  general  integer  variables;  exact  algorithm; 
ceiling  points;  implicit  enumeration;  linear  programming  relaxation 
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This  report  describes  an  exact  algorithm  for  solving  the  pure,  general  integer  linear 
programming  problem  in  m  constraints  and  n  variables  xj,  j  =  1,  whose  form  is 

Maximize  cTx  =  z 

subject  to  Ax  <b  ( ILP ) 

x  >  0,  x  integer, 


where  A  G  b  €  §im  and  c  G  X”.  All  the  data  {A,6,c}  are  assumed  to  be  rational 

numbers,  but  they  are  unrestricted  in  sign.  The  problem  is  “pure”  in  that  all  of  the 
variables  are  required  to  take  on  nonnegative  integer  values.  It  is  “general”  in  the  sense 
that  the  variables  may  take  on  any  nonnegative  integer  values  permitted  by  Ax  <  b,  as 
opposed  to  being  restricted  to  0  or  1  (the  binary  case).  An  important  additional  assumption 
is  that  no  implicit  or  explicit  equality  constraints  are  used  to  define  the  feasible  region 
FR  =  {x  >  0|  Ax  <  6}  for  ( LPr ),  the  linear  programming  relaxation  associated  with 
(ILP).  Common  applications  of  this  model  occur  in  capital  budgeting  (project  selection), 
resource  allocation  and  fixed-charge  (plant  location)  problems.  A  further  discussion  of 
application  areas  for  (ILP)  may  be  found  in  [11]  or  [26]. 

Another  report  by  the  authors  [24]  describes  a  procedure  for  approximately  solving 
(ILP)  called  the  Heuristic  Ceiling  Point  Algorithm.  Like  all  heuristic  algorithms,  it  offers 
no  guarantee  of  finding  an  optimal  solution.  In  fact,  while  the  Heuristic  Ceiling  Point 
Algorithm  usually  finds  solutions  of  very  hiqh  quality,  the  algorithm  is  not  even  guaran¬ 
teed  to  find  a  feasible  solution.  Furthermore,  even  with  an  optimal  integer  solution  at 
hand,  the  Heuristic  Ceiling  Point  Algorithm  can  only  rarely  verify  (prove)  that  such  a 
solution  is  indeed  optimal.  For  these  reasons,  an  exact  algorithm  was  developed  based  on 
a  more  systematic  search  for  feasible  “1-ceiling  points.”  To  understand  why  this  might  be 


a  reasonable  approach,  we  briefly  review  some  of  the  key  concepts  of  [23]. 

An  integer  solution  x  is  a  1-ceiling  point  with  respect  to  the  ith  constraint,  denoted 


x  =  l-CP(i),  if  (1)  x  satisfies  this  constraint,  re.,  ofx  <  bi  (where  a*  is  the  ith  row  of - - 

_ i 

the  constraint  matrix  A),  and  (2)  modifying  some  component  of  x  by  +1  or  —1  yields  a  ’odes  j 


□  □ 
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solution  which  violates  this  constraint,  i.e.,  ajx  +  1  >  6,  for  at  least  one  j.  Thus,  x  = 

l-CP(t)  means  x  narrowly  satisfies  the  ith  constraint:  taking  a  unit  step  from  x  toward 
the  ith  constraining  hyperplane  in  a  direction  parallel  to  some  coordinate  axis  results  in 
an  infeasible  point.  Similarly,  an  integer  solution  x  is  defined  to  be  a  1-ceiling  point  with 
respect  to  the  feasible  region  FR,  denoted  x  =  l-C  P(FR),  if  (1)  x  satisfies  all  constraints 
and  (2)  modifying  some  component  of  x  by  -fl  or  —  1  leads  to  a  solution  which  violates 
one  or  more  constraints,  i.e.,  3i  :  aTx  +  lay  |  >  bi  for  at  least  one  j.  It  is  demonstrated 
in  [23]  that  all  optimal  solutions  for  an  ( ILP )  whose  feasible  region  is  non-empty  and 
bounded  are  feasible  l-CP(i)’s,  i.e.,  l-C  P(FR)'s.  Consequently,  one  way  to  solve  (ILP) 
is  to  enumerate  its  feasible  1-ceiling  points.  Our  heuristic  approach  is  based  upon  the  idea 
that  a  feasible  1-ceiling  point  found  relatively  near  x  is  apt  to  have  a  high  (possibly  even 
optimal)  objective  function  value.  On  48  test  problems  taken  from  the  literature,  searching 
for  such  1-ceiling  points  usually  did  provide  a  very  good  solution  with  a  moderate  amount 
of  computational  effort. 

The  Exact  Ceiling  Point  Algorithm  begins  by  executing  the  Heuristic  Ceiling  Point 
Algorithm  (which  includes  solving  the  linear  programming  relaxation),  so  Section  2  ex¬ 
amines  the  assumptions  needed  for  these  steps.  The  goal  at  this  stage  is  to  be  able  to 
construct  a  small  but  sufficient  search  region.  Section  3  discusses  how  this  search  region 
is  divided  up  to  permit  the  search  to  focus  on  finding  1-ceiling  points  with  respect  to  one 
specific  constraint.  This  leads  to  the  calculation  of  unconditional  variable  bounds  which 
correspond  geometrically  to  a  search  hyperrectangle.  The  overall  process  of  searching  for 
1-ceiling  points  within  this  search  hyperrectangle  is  outlined  at  the  beginning  of  Section  4 
and  then  each  of  three  subsections  provide  more  detail  on  a  particular  aspect  of  the  search. 
The  first  subsection  describes  the  calculation  of  conditional  variable  bounds  which  are  at 
least  as  strong  as  the  unconditional  variable  bounds;  the  second  subsection  describes  an 
efficient  aspect  of  the  enumeration  procedure  called  “double  backtracking”  which  allows 
subregions  of  the  search  rectangle  to  be  skipped  altogether;  the  third  subsection  describes 
what  occurs  when  a  new  incumbent  solution  is  found.  An  overview  of  the  entire  Exact 
Ceiling  Point  Algorithm  is  given  is  Section  5.  Section  6  describes  a  preliminary  search 
procedure  on  a  constraint  called  an  “intersection  cut,”  where  this  procedure  is  designed 
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to  speed  up  our  algorithm  by  searching  the  region  closest  to  x  first.  The  choice  of  a 
search  constraint  in  the  Exact  Ceiling  Point  Algorithm  is  discussed  in  Section  7,  while 
a  few  other  computational  issues  are  examined  in  Section  8.  Section  9  reports  on  our 
computational  experience.  Section  10  offers  some  conclusions  about  our  approach  and  is 
followed  by  two  appendices.  The  first  appendix  gives  the  variable  bounds  and  options  used 
in  the  GAMS/ZOOM  runs  reported  in  Section  9,  while  the  second  lists  the  Fortran  code 
implementation  of  the  Exact  Ceiling  Point  Algorithm. 


2.  Assumptions 


The  Exact  Ceiling  Point  Algorithm  begins  by  finding  an  optimal  solution  x  for  the  LP- 
relaxation  ( LPr ),  the  set  A  of  constraints  binding  at  x,  and  the  set  of  extreme  directions 
emanating  from  x  which  form  the  cone  FR.  The  Exact  Algorithm  also  uses  an  initial 
feasible  integer  solution,  xH,  provided  by  the  Heuristic  Algorithm.  The  assumptions  made 
here  are  very  similar  to  those  specified  in  [16]. 

Assumption  1:  The  set  of  feasible  solutions  for  (LPr)  is  non-empty  and  bounded. 
Assumption  2:  The  optimal  solution  found  for  (LPr)  is  not  all-integer. 

Assumption  3:  The  unique  optimal  solution  for  (LPr)  is  x. 

Assumption  4:  A  feasible  solution  xr  for  (ILP),  with  objective  value  z_,  is  known. 

The  first  assumption  implies  that  f  exists.  The  second  implies  that  (ILP)  is  not 
solved  simply  by  solving  (LPr),  so  that  we  have  a  need  for  an  exact  algorithm.  The 
third  assumption  is  the  most  serious;  however,  when  x  is  not  unique,  there  are  ways  to 
perturb  the  data  of  (ILP)  without  altering  its  optimal  solution(s)  so  that  this  condition 
does  hold.  Another  alternative  is  to  continue  to  append  cutting  planes  to  the  problem 
until  a  unique  x  is  found  (see  [16,  pp.  670-673]).  The  last  assumption  is  not  really  needed 
for  the  Exact  Ceiling  Point  Algorithm  to  run,  but  vastly  improves  its  performance  since 
it  provides  another  relatively  strong  constraint,  cTx  >  z,  on  the  solutions  that  need  to 
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be  considered.  If  no  feasible  integer  solution  for  ( ILP )  has  been  identified  prior  to  the 
start  of  the  Exact  Algorithm,  then  the  algorithm  starts  with  z  =  —  oo.  Alternatively,  the 
value  of  some  feasible  non-integer  solution  could  be  used  sis  a  lower  bound  on  the  optimal 
objective  value  of  (ILP).  Then,  if  it  is  found  that  there  do  not  exist  any  feasible  integer 
solutions  whose  objective  value  exceeds  or  equals  this  bound,  so  the  bound  is  not  a  valid 
one,  the  Exact  Algorithm  would  have  to  be  restarted  with  a  smaller  lower  bound,  and  the 
previous  lower  bound  would  then  become  an  upper  bound  [16,  p.  673]. 

With  one  more  definition,  we  will  be  able  to  define  a  volume  which  must  be  searched 
to  find  an  optimal  solution  for  (ILP).  For  k  =  l,...,n,  let  pk  be  defined  as  the  point 
of  intersection  between  the  kth  extreme  direction  dk  emanating  from  i  and  the  objective 
constraint  hyperplane  cTx  =  z.  These  intersection  points  {pfc,  k  =  1, ...,  n}  exist  as  long  as 
x  exists  and  the  objective  constraint  hyperplane  is  not  parallel  to  any  edge  or  face  of  the 
cone  FR  (Assumptions  1  and  3).  The  point  pk  has  the  form  x  +  A kdk,  where  A*  is  such 
that 

cT(x  +  A  kdk)  =  z. 


Solving  for  A*  yields 


A*  =  (z  -  cTx)/cTdk. 


Therefore, 


pk  =  x  +  (z  —  cTx)dk/crdk,  for  k  =  l,...,n. 


For  notational  convenience  in  what  follows,  define  p°  =  x.  Now  we  can  restate  Hillier’s 
key  theorem  [16,  Theorem  1]  for  this  point  in  the  analysis. 


Theorem  1.  Under  Assumptions  1,  2,  and  3,  all  optimal  solutions  for  (ILP)  are  contained 
in  the  n-simplex  S  whose  (n  +  1)  extreme  points  are  {p° ,pl,  ...,pn}. 

An  example  in  3?3  of  an  n-simplex  5  is  shown  in  Figure  1,  where  the  extreme  directions 
emanating  from  i  are  {d1 ,  d2,  d3}.  Theorem  1  indicates  that  once  z_  has  been  specified, 
our  attention  may  be  confined  to  S  without  fear  of  missing  any  optimal  solutions  for 
(ILP).  Consequently,  any  non-binding  constraint  which  does  not  intersect  5  may  be 
disregarded  henceforth.  In  the  Exact  Algorithm,  any  constraint  which  is  not  violated  by 
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Figure  1.  An  n-simplex  S  in  J?3  with  vertices  {x.p^p^p3}  and  the  corresponding  search 
hyperrectangle  SHR(S). 

at  least  one  pk  is  redundant  and  dropped  from  the  further  consideration.  On  the  other 
hand,  the  n-simplex  may  be  larger  than  necessary  due  to  the  presence  of  (non- redundant) 
non-binding  constraints  which  chop  off  parts  of  the  n-simplex.  Notice  that  if  we  seek 
only  solutions  which  are  strictly  better  than  x#  and  all  the  cj  are  integer,  the  objective 
constraint  hyperplane  to  use  in  the  determination  of  {pfc,  k  =  l,...,n}  is  cTx  =  :  +  1.  It 
seems  worthwhile  to  recompute  this  set  of  intersection  points  every  time  a  new  (higher¬ 
valued)  incumbent  solution  is  found  since  this  reduces  the  size  of  the  n-simplex.  Finally, 
we  may  calculate  a  set  of  lower  and  upper  bounds  [5^,5^]  for  Xj  by  finding  the  minimum 
and  maximum,  respectively,  over  the  jth  component  of  all  vertices  of  the  n-simplex  and 
rounding  appropriately.  The  lower  bound  S_j  is  not  permitted  to  be  less  than  0: 


Sj  =  max{fmin{pf}],0},  for  j  =  l,...,n, 
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and 

=  |max{pJ}J ,  for  j  =  1, . . . , n. 

These  bounds  define  what  we  shall  refer  to  as  SHR{S ),  the  search  hyperrect angle  for  the 
n-simplex  S ,  as  illustrated  in  Figure  1.  If,  for  any  component  j,  we  find  that  the  upper 
bound  5j  is  strictly  less  than  the  the  lower  bound  §_j,  then  SHR(S)  —  0.  In  this  case,  the 
problem  is  solved  because  there  are  no  feasible  integer  solutions  better  than  the  incumbent. 
Otherwise,  we  begin  to  search  for  l-CP(FP)’s. 


3.  Unconditional  Variable  Bounds 


The  preceding  section  indicated  that  to  solve  (JLP)  it  is  sufficient  to  enumerate  the 
feasible  1-ceiling  points  contained  in  the  region  defined  by  SHR(S).  To  improve  effi¬ 
ciency,  we  split  up  the  n-simplex  and  further  confine  our  search  to  subhyperrectangles  of 
SHR(S).  Like  the  Heuristic  Ceiling  Point  Algorithm,  the  Exact  Algorithm  seeks  1-ceiling 
points  with  respect  to  one  search  constraint  at  a  time.  Therefore,  a  search  hyperrectangle 
SHR(i )  for  constraint  (t)  contained  within  SHR(S)  is  constructed  which  is  large  enough 
to  contain  all  feasible  l-CP(i)’s  with  objective  value  greater  than  z.  This  subhyperrect¬ 
angle  is  constructed  once  a  promising  search  constraint  has  been  identified,  which  is  the 
topic  of  Section  7. 

FVom  Definition  3.4,  a  necessary  and  sufficient  condition  for  an  integer  solution  x  to 
be  a  1-ceiling  point  with  respect  to  constraint  (i)  is 

x  =  1— CP(t)  «=>  0  <  ft,-  -  ajx  <  max  |ajj|  -  1  (1) 

where  all  coefficients  aij  are  assumed  to  be  integer.  (If  not  all  coefficients  fire  integer, 
subtract  a  small  positive  quantity  e  instead  of  1  from  the  right  hand  side  of  (1).  This  will 
be  the  case  when  we  search  the  intersection  cut  described  in  Section  6.)  In  other  words, 
for  i  to  be  a  1-ceiling  point  with  respect  to  constraint  (*),  it  must  satisfy  the  constraint 


ofx  <  bi 


(*) 
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but  not  be  too  far  away  from  the  i<A  constraint  hyperplane,  t.  e. ,  satisfy 

o[x>bi-ti  (*') 

where  t{  =  max;  |a^j  —  1  may  be  thought  of  as  the  distance  constraint  (*)  is  translated  in 
order  to  impose  the  ceiling  point  condition.  Thus,  only  those  integer  solutions  satisfying 
both  constraints  (i)  and  ( i ')  are  l-CP(i)’s.  These  two  constraints  are  central  to  the 
construction  of  the  appropriate  search  region.  An  example  in  3ft2  is  given  in  Figure  2. 


Figure  2.  All  feasible  l-CP(l)’s  with  objective  value  as  large  as  £  lie  in  the  part  of  S 
between  (1)  and  (1'). 


Now  let  Ei  be  the  set  of  extreme  rays  emanating  from  x  which  lie  on  the  ith  constraint 
hyperplane.  Also  let  qk  be  the  point  of  intersection  between  the  kth  extreme  ray  of  FR 
and  the  translated  constraint  hyperplane  (i'),  for  all  k  Ei.  In  Figure  2,  the  point  q2  is 
shown  as  the  intersection  of  constraint  hyperplane  (1')  and  extreme  ray  2,  where  extreme 
ray  2  coincides  with  constraint  hyperplane  (2).  When  exactly  n  constraints  axe  binding  at 
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x,  each  translated  constraint  hyperplane  (t')  intersects  just  one  extreme  ray  because  the 
other  (n  —  1)  extreme  rays  he  on  constraint  hyperplane  (t)  which  is  parallel  to  (t').  In  the 
case  where  *  is  overdetermined,  more  than  one  extreme  ray  will  intersect  the  translated 
constraint  hyperplane  (»').  The  point  g*  has  the  form  x  +  7 kdk,  where  7*  is  such  that 

of(x  +  ikdk)=:bi-ti. 


Solving  for  7*  yields 

7*  =  ~ti/ajdk. 

Therefore, 

qk  =  x-(ti/aJdk)dk,Vk<tEi. 

The  search  hyperrectangle  for  this  constraint,  SHR(i ),  is  defined  by  the  ranges 
{[/j,  Uj],  j  =  1, ..., n},  where  the  bounds  for  each  x}  are  formed  by  taking  the  minimum  and 
maximum,  respectively,  over  the  jth  component  of  all  of  the  intersection  points  required 
to  define  the  hyperrectangle’s  vertices  and  rounding  appropriately.  Furthermore,  these 
bounds  should  be  no  wider  than  those  defining  SHR(S)  because  searching  beyond  the 
boundary  of  SHR(S)  is  unproductive: 

lj  =min{[min{min{pj},  min{g*}}],  S} }  (2/) 

K%C*% 

and 

Uj  =  max{|max{max{p{},  mf«{g*}}J,  ?,}.  (2u) 

If,  for  any  component  j,  we  find  uj  <  lj,  then  SHR(i)  =  0,  implying  that  there  are 
no  integer  solutions  at  all  in  this  search  hyperrectangle,  and  a  new  search  constraint  is 
identified.  Assuming  this  is  not  the  case,  we  proceed  to  enumerate  the  integer  solutions 
contained  within  the  search  hyperrectangle,  as  described  in  the  next  section. 


4  Enumerating  Solutions  Within  a  Search  Hyperrectangle _ £ 

4.  Enumerating  Solutions  Within  a  Search  Hyperrectangle 

Once  the  search  hyperrectangle  SHR(i )  has  been  defined  for  search  constraint  (t), 
finding  l-C'P(t)’s  better  than  the  incumbent  amounts  to  examining  SHR(i)  for  solutions 
which  are  feasible  with  respect  to  all  the  relevant  constraints,  including  ( i ').  Because 
this  search  hyperrectangle  is  contained  within  the  n-simplex,  i.e.,  SHR(i)  C  SHR(S), 
any  feasible  integer  solution  found  will  become  the  new  incumbent  solution.  The  overall 
enumeration  process,  shown  as  a  flow  diagram  in  Figure  3,  contains  three  features  which 
differentiate  it  from  a  simple  exhaustive  enumeration  scheme:  the  use  of  conditional  vari¬ 
able  bounds,  a  potential  for  “double  backtracking,”  and  a  tightening  of  bounds  when  a 
new  incumbent  is  found.  These  are  described  in  each  of  the  next  three  subsections. 

4.1.  Conditional  Variable  Bounds 

Suppose  the  variables  are  rearranged  so  that  Xj  is  fixed  first,  X2  is  fixed  second,  and  so 
forth.  Given  a  “partial  solution”  (xj,x2,...,xJ_i)  whose  components  are  fixed  in  value,  a 
“completion”  of  the  partial  solution  refers  to  an  assignment  of  values  to  the  remaining  free 
components.  Fixing  the  value  of  the  first  variable  Xj  to  an  integer  value  between  Zj  and  ux 
may  reduce  the  range  of  values  for  some  or  all  of  the  free  variables  (x2,  ...,x„).  Similarly, 
fixing  the  first  ( j  —  1)  variables  may  tighten  the  bounds  on  the  remaining  variables.  In 
other  words,  it  is  possible  to  calculate  conditional  bounds  on  the  free  variables  given  the 
value  of  the  fixed  variables. 

As  in  the  Heuristic  Ceiling  Point  Algorithm,  we  assume  that  all  constraints  are  in 
<  form,  having  multiplied  any  >  constraints  through  by  —1  if  necessary.  The  following 
procedure  develops  conditional  bounds  [Lj\xj-i,Uj\xj-x]  for  the  variable  Xj,  given  that 
variables  xj,  ...,Xj_j  are  fixed  in  value  and  x^, ...,  x„  are  free.  (In  Figure  3  these  conditional 
bounds  are  abbreviated  as  [Lj,Uj].)  First,  some  useful  notation  is  introduced.  Let 

j- 1 

9ij  =  b,  —  y  ^  OijX k 

k=  1 
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^  Initializations^ 

1 

t 

Next  va 
X]  «  x 

lue  of  Xj 

1  ♦  axj 

Forward  5tep 
J- J«  1 


t 


Calculate/update 
conditional  Pounds 
(LJ.  UJ]  on  XJ  (I) 


i  c  2™P 


Figure  3.  Flow  diagram  of  the  Exact  Algorithm's  enumeration  scheme. 

(1)  See  Section  4.1  for  further  explanation. 

(2)  See  Section  4.2. 

(3)  See  Section  4.3. 


4  Enumerating  Solutions  Within  a  Search  Hyperrectangle  _  II 


be  the  gap  or  alack  remaining  in  constraint  (i)  after  having  fixed  x% , ...,  Xj-\.  Also,  let 

n 

mm{aiklk,aikuk}  (3) 

k=j+l 

be  the  minimum  amount  of  the  gap  that  may  be  used  up  by  fixing  variables  xJ+1 , xn 
within  their  unconditional  bounds.  Feasibility  of  a  complete  solution  x  with  respect  to  the 
ith  constraint  requires 


j-l  n 

Y2aikXk  +  ciijXj  +  ^ 

aikXk  <  bi 

fc=i  *==j+i 

i-i 

n 

CLijXj  <  bi  ^  ^  aikxk 

y.  aikXk 

i=i 

*=>+ 1 

n 

aijXj  <  9ij  -  ^2 

aikXk 

t=i+i 

=>  “a*]  <  9ij  -  (4) 

If  a,;  is  positive,  then  (4)  yields  an  upper  conditional  bound  for  x}  when  just  constraint 
(i)  is  considered: 

x,  <  (gij  -WiJ+l)/aij.  (5) 

On  the  other  hand,  if  CLij  is  negative,  (4)  yields  a  lower  conditional  bound  for  Xj  when  just 
constraint  (i)  is  considered: 

xj  >(9ij  ~witj+i)/aij.  (6) 

Note  that  the  numerator  may  be  negative,  so  that  the  lower  conditional  bound  is  not 
necessarily  negative.  If  a,j  is  positive  for  all  (t),  (4)  yields  only  upper  conditional  bounds. 
Finally,  if  a,,  =  0,  constraint  (:)  cannot  be  used  to  determine  a  conditional  bound  on  x}. 
Letting  Lij\xj-X  and  denote  the  lower  and  upper  conditional  bounds,  respectively, 

on  Xj  when  constraint  ( i )  is  considered  alone  after  fixing  variables  (xj, ...,  Xj_i),  we  have 


T  _  /  (ffij  Wi,j+i  )/aiji  if  a*j  <  0> 

X'ijlXj-i  =  | ifao  >0, 

Cf  lx  ,  =/“>’  if  <  °; 


and 
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where  lj  and  uj  are  the  unconditional  lower  and  upper  bounds  on  xy  defined  earlier.  To 
maintain  feasibility  with  respect  to  all  constraints,  the  conditional  lower  and  upper  bounds 
\Lj\xj-i >Uj\xj-i)  on  Xj  must  be  the  tightest  among  all  those  conditional  bounds  which 
consider  only  one  constraint  at  a  time.  Therefore,  the  desired  conditional  lower  and  upper 
bounds  are 


Li\xi- 1  =  rmax{I0|i>_1}]  (7) 

and 

Uj\xi- 1  =  Lmin{f/,>|ary_, }  J .  (8) 

When  we  find  Lyjxy_i  >  C/y|xy_i,  there  are  no  feasible  completions  of  the  partial  solution 
(*i,  —, In  this  case,  we  backtrack  to  the  next  value  of  xy_j.  The  bounds  (7)  and  (8) 
are  vey  similar  to  those  defined  in  a  lemma  by  P.  Krolak  which  are  central  to  his  Bounded 
Variable  Algorithm  [18]. 

At  first  glance,  it  may  appear  that  calculating  £,,(*,_!  and  Cf,y|xy_!  each  time  xy_j 
changes  is  computationally  prohibitive,  but  actually  it  is  not.  This  is  because  L^xj^ 
and  Uij\xj-i  both  change  in  a  predictable  way  as  xy_i  changes.  In  fact,  each  is  a  linear 
function  with  a  slope  that  can  be  calculated  once  at  the  outset  of  the  Exact  Ceiling  Point 
Algorithm,  as  we  now  demonstrate.  To  see  the  effect  of  changing  xy_ lt  let  vj-i  be  its 
current  value  and  v'_j  the  value  to  which  it  is  then  changed,  so  t>'_,  =  vy_ j  +  6j-i, 
where  tfy.j  €  {-1,+1}.  (Our  rule  is  to  set  =  -1  if  c,_,  >  0,  or  set  6j-\  =  +1  if 
Cj-\  <  0,  implying  that  we  start  at  the  more  attractive  end  of  the  interval  defined  by  the 
unconditional  botrnds.)  Let  Lij\x'j_1  and  Uij\x'j_1  denote  the  lower  and  upper  conditional 
bounds,  respectively,  on  ij  when  considering  constraint  (i)  alone  after  fixing  variables 
(x1,...,Xj_2,xJi_j)  to  (oi,...,t;y_2,»y_I).  Similarly,  let  ffij\xj-i  denote  the  quantity  gij 
given  the  partial  solution  (t/j , ...,  Vj-2,  ry_i )  and  y,y|x'_j  the  quantity  gij  given  the  partial 
solution  (»i,..., vy_2,v'-_ j).  Assuming  that  o^-  >  0,  we  are  interested  in  how  Uij\xj-\ 
differs  from  Uij\x'j_x: 
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U*j\xj-i  Uij\xj-i  —  Ksr«l*i-1  ~  w*j+i)  —  ~  u,«,i+i)]/°«i 

=  (9ij\xj-l  ~  9ij\xj—l)/aij 
=  (-ai,i-iv>-i  -  (-a»,i-iV>-i)]Mi 
=  ai,j-i(vi-i  ~  v'j-i)/aij 

where  fij  =  Since  both  fij  and  6j~\  are  independent  of  the  value  of  Xj_i, 

^«ilxj-i  is  a  linear  function  of  Xj-\  with  slope  ifij.  When  a^  <  0,  a  similar  derivation 
reveals  that  Lij\xj-\  is  also  a  linear  function  of  x;_j  with  slope  ±/tJ.  Therefore,  after 
calculating  U%j\xj-\  and  Lij\xj-\  for  the  initial  value  of  x^_i  within  its  conditional  bounds, 
one  of  the  two  conditional  bounds  always  changes  in  an  additive  fashion  while  the  other 
remains  equal  to  the  value  of  its  unconditional  bound: 


and 


■^olxj-i  = 

KM- 1  = 


w 


{  ui, 

\  —  fij  —  1  fij » 


if  aij  <  0; 
if  aij  >  0, 


if  aij  <  0; 
if  aij  >  0. 


4.2.  Double  Backtracking 


In  addition  to  being  easy  to  update,  the  quantities  Lij\xj-i  and  |xy_i  impart  an 
important  property  to  the  conditional  bounds  Lj\xj-\  and  Uj\xj-X .  As  the  maximum  over 
a  set  of  Unear  functions,  Lj\x^x  is  a  (piecewise  linear)  convex  function  of  x;_i  by  a  well- 
known  theorem  (see  [3,  Theorem  4.13]  for  example).  Similarly,  Uj\x,-i  is  the  minimum 
over  a  set  of  linear  functions  and  therefore  is  a  (piecewise  linear)  concave  function  of  x3  _ , . 
This  is  illustrated  in  Figure  4  for  the  case  in  which  ai}  >  0,  for  all  (i,j),  so  that  the 
conditional  lower  bound  Lj\xj-t  is  equal  to  the  unconditional  boimd  lj. 

As  ij_ i  ranges  between  its  conditional  bounds  Lj-i  |x/_2  and  Uj-i  |x;_2,  the  condi¬ 
tional  bounds  on  x;,  L;|xj_i  and  l/j|xy_j ,  may  cross  once,  twice  or  not  at  all.  If  we  observe 
that  Lj\ij-i  exceeds  Uj\xj-\  for  a  particular  value  of  xj- j,  we  can  certainly  backtrack  to 


i 


Figure  4.  Upper  and  lower  conditional  bounds  on  Xj  when  all  ai7  >  0. 

the  next  value  of  Xj-\ ,  because  there  are  no  feasible  completions  of  the  partial  solu¬ 
tion  (*|,. . .  ,Xj-i).  But  there  is  also  a  possibility  that  wifi  exceed  Uj\xj-\  for  all 

remaining  values  of  Xj_i,  due  to  the  convexity  and  concavity  properties  of  the  conditional 
bounds.  In  this  case,  we  can  “double  backtrack,”  *.e.,  backtrack  to  the  next  value  of  Xj_2. 
The  conditions  for  when  it  is  possible  to  double  backtrack,  thereby  significantly  reducing 
the  amount  of  enumeration,  are  given  in  the  next  theorem. 

Theorem  2.  Suppose  the  backtracking  condition  Lj\xj-\  >  Uj\xj-i  holds  for  a  particu¬ 
lar  value  of  Xj- 1,  implying  that  there  are  no  feasible  completions  of  the  partial  solution 
(xi,. . .  ,Xj_i).  A  sufficient  condition  to  double  backtrack,  *.e.,  skip  over  all  the  remaining 
integer  values  of  Xj-\  within  its  conditional  bounds,  is  that 

for  some  iL  g  argmax;{L;;jx,_i }  and  some  iv  g  argmin,{C/tJ|x^_i}. 
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Proof: 

Lj\x'i- x  >  Lj\xj- 1  -  fiLjSj-u  where  iL  £  argmax^ix^j} 

>  Uj\ij-i  —  fiLjSj~ |,  by  the  backtracking  condition 

>  Uj\xj-i  -ficrj5j- 1,  where  ty  €  argn^nf^lx^i} 

> 

The  first  and  last  inequalities  hold  by  the  convexity  and  concavity  of  L} |x'_x  and  Uj |x'_x, 
respectively,  while  the  third  inequality  follows  if  the  double  backtracking  condition  stated 
by  the  theorem  holds.  Repeating  the  argument  for  subsequent  values  of  Xj-\ ,  we  find  that 
the  lower  conditional  bound  on  x;  exceeds  the  corresponding  upper  conditional  bound  not 
only  for  x' _x  but  for  all  remaining  values  of  x'_x  within  the  conditional  bounds  of  xy_x. 
Thus,  we  can  double  backtrack.  I 

Condition  (9)  is  stronger  than  necessary  because  there  may  be  some  other  constraint 
(i')  which  is  not  a  member  of  the  argmaXj{L,^|xj_i }  but  which  is  a  member  of  the 
argmax,-{£tj|x'-_1}  so  that  Lj\x'j_l  >  Lj\xj-i  —  fiijtj- 1-  However,  constraint  indices 
il  and  iv  are  readily  available,  having  been  identified  in  the  calculation  of  Lj\xj-i  and 
respectively.  Note  that  it  only  makes  sense  to  check  the  double  backtracking  con¬ 
dition  after  we  have  found  that  we  can  perform  an  ordinary  backtracking  step.  We  must 
check  the  condition  of  Theorem  2  because  of  the  possibility  that  the  conditional  bounds 
may  cross  twice  as  x;_x  moves  between  its  conditional  bounds.  If  we  double  backtracked 
after  the  first  observation  of  the  lower  conditional  bound  exceeding  the  upper  conditional 
bound,  we  could  miss  some  partial  solutions  with  feasible  completions.  Double  backtrack¬ 
ing  seems  to  be  most  helpful  on  problems  in  which  the  unconditional  bounds  are  fairly 
wide  for  one  or  more  variables,  such  as  the  fixed-charge  problems  discussed  in  Section  9. 

4.3.  When  a  New  Incumbent  is  Found 

The  third  feature  of  our  Exact  Algorithm  which  distinguishes  it  from  ordinary  enu¬ 
meration  is  manifest  when  a  feasible  solution  is  completed.  Because  all  constraints  axe 
taken  into  account  in  the  calculation  of  the  conditional  variable  bounds,  including  an 
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objective  function  constraint,  any  completed  solution  is  not  only  feasible  but  also  has 
a  strictly  superior  objective  function  value  to  that  of  the  current  incumbent.  As  such, 
it  becomes  the  new  incumbent  solution.  At  this  stage,  the  intersection  points  between 
the  extreme  rays  emanating  from  x  and  the  new  objective  function  constraint  hyperplane 
(which  define  the  n-simplex)  move  closer  to  x.  Thus,  the  n-simplex  shrinks  and  a  tighter 
set  of  integer  component  bounds  defining  the  SHR(S)  may  result.  If  so,  the  unconditional 
variable  bounds  developed  for  the  current  search  constraint  (»)  will  also  tighten,  that  is, 
move  closer  together.  Consequently,  the  minimal  completion  values  tv,j  defined  in  (3)  will 
change,  possibly  leading  to  tighter  conditional  variable  bounds.  Thus,  as  we  backtrack 
from  the  completed  solution  to  new  partial  solutions,  new  conditional  variable  bounds  are 
computed  which  reflect  the  improved  objective  value  of  the  incumbent.  The  magnitude  of 
the  change  in  the  objective  function  value  determines  the  extent  to  which  these  modified 
conditional  bounds  speed  up  the  enumeration  process. 


5.  Overview  of  the  Exact  Ceiling  Point  Algorithm 


Having  seen  how  to  enumerate  1-ceiling  points  with  respect  to  a  specific  constraint, 
we  are  now  in  a  position  to  see  where  this  process  fits  into  the  overall  Exact  Ceiling  Point 
Algorithm  described  in  Figure  5.  The  iterative  part  (Step  3)  of  the  Exact  Algorithm 
essentially  “branches”  on  one  of  the  functional  constraints  (t)  and  seeks  to  “fathom”  a 
subregion  SHR(i )  of  the  SHR(S )  by  enumerating  implicitly  or  explicitly  all  l-CP(i)’s 
within  SHR(i).  Suppose  that  the  best  feasible  integer  solution  known  after  having  searched 
this  region  is  Xi  whose  objective  function  value  is  2  =  cTX{.  If  the  new  n-simplex  S'  is 
formed  by  the  objective  function  hyperplane  cTx  =  2  +  1  and  the  constraints  binding  at 
x,  then  S'  excludes  the  best  known  solution  x<.  It  is  possible,  then,  that  S'  contains  no 
integer  solutions  at  all.  A  sufficient  condition  for  this  is  if  there  is  some  j  such  that  the 
new  upper  bound  3^  is  strictly  less  than  the  new  lower  bound  §_j.  In  this  case,  the  problem 
has  been  solved.  If  not,  we  replace  the  previous  search  constraint  (i)  by  a  constraint  ( i ") 
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parallel  to  (i)  whose  right  hand  side  is  —  maxj  |ajj-|,  resulting  in  a  modified  problem 

(ILP1)  which  remains  to  be  searched.  The  modified  relaxation,  ( LP'R ),  is  then  solved  for 
its  optimal  solution  and  value  ( x\z ').  Note  that  z'  <  z  and  that  z'  is  an  upper  bound  for 
the  optimal  objective  function  value  for  ( ILP ').  Therefore,  if  z'  <  z,  the  problem  is  solved 
with  x ,  as  an  optimal  solution  for  (ILP).  If  the  problem  is  not  solved,  a  new  iteration 
begins  by  selecting  another  search  constraint  and  enumerating  its  1-ceiling  points.  Thus, 
the  algorithm  proceeds  in  a  “search  and  cut”  fashion  until  the  problem  is  solved.  (Step  2 
will  be  described  in  the  next  section.) 

Step 

0.  a.  Solve  ( LPr )  =>  (x,z  =  cTx ). 

b.  Apply  Heuristic  Ceiling  Point  Algorithm  =►  (xH,z_~  cTxH). 

1.  Construct  n-simplex  S  and  associated  search  hyperrectangle  SHR(S). 

2.  a.  Construct  Intersection  Cut  (I)  and  search  hyperrectangle  SHR(I). 

b.  Enumerate  feasible  l-CP(J)’s  within  SHR(I).  =>  (i/,  r  =  cTi/). 

c.  Reduce  J  to  max* {cTrk},  where  r*  s  intersection  of  kth  extreme  ray  and  (/). 

d.  Stop  if  z  =  7;  otherwise, 

3.  REPEAT 

a.  Select  a  search  constraint  (i)  and  construct  SHR(i). 

b.  Systematically  enumerate  feasible  l-CP(i)’s  =>  (x,,z  =  crx,). 

c.  Replace  of x  <  bi  with  afx  <bi  ~  maxj  |a,;|. 

d.  Resolve  ( LP'R )  to  obtain  (x\  z*  =  cTx'). 

UNTIL  OPTIMAL  (FR  =  0  or  z  >  z'). 


Figure  5.  Outline  of  the  Exact  Ceiling  Point  Algorithm. 

Theorem  3.  Suppose  the  set  of  feasible  solutions  for  (ILP)  is  non-empty  and  bounded 
(Assumption  1),  x  is  unique  (Assumption  3)  and  z  >  —  oo.  Then  the  iterative  procedure 
(Step  3)  of  the  Exact  Ceiling  Point  Algorithm  given  in  Figure  5  is  guaranteed  to  find  an 
optimal  solution  for  (ILP)  in  a  finite  number  of  steps. 
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Proofs  By  Theorem  2  of  [23]  and  Theorem  1  above,  the  total  region  T  which  needs  to  be 
examined  in  order  to  solve  the  problem  is  the  portion  of  the  n-simplex  S  which  contains 
all  of  its  l-CP(FR)'s: 

T  =  |J{*  >  0|  bi  -  max  |o0|  <  afx  <  bi,  cTx  >  z}. 

i 

Since  T  consists  of  subsets  of  S,  T  itself  is  bounded.  There  are  only  a  finite  number  of 
constraints,  each  of  which  is  searched  at  most  once  for  its  1-ceiling  points.  What  remains 
to  be  shown  is  that  the  scheme  of  Section  4  for  enumerating  l-CP(»)’s  with  respect  to  a 
specific  constraint  (i)  is  finite. 

For  a  given  constraint  (i),  a  search  hyperrectangle  SHR(i )  is  defined  and  then  exam¬ 
ined  for  l-CP(t)’s.  Being  a  subset  of  T,  this  search  hyperrectangle  $HR(i)  is  bounded, 
implying  that  all  of  the  corresponding  unconditional  variable  bounds  j  =  1, . . .  ,n} 

are  bounded.  The  conditional  variable  bounds  {[Lj\xj-t,  Uj\xj^t],  j  =  l,...,n}  are  also 
bounded  since,  by  definition,  >  lj  and  lfj\xj-i  <  Uj  for  all  j.  Thus,  while  the  total 

number  o,'  integer  solutions  contained  within  SHR(i)  may  be  large,  it  is  finite.  Finally, 
the  scheme  which  enumerates  integer  solutions  within  SHR(i)  examines  at  most  once  each 
complete  solution  or  partial  solution  leading  to  one  or  more  complete  solutions.  Partial 
solutions  are  fathomed  if  they  are  shown  to  be  unable  to  lead  to  a  feasible  completion  that 
is  better  than  the  incumbent.  Any  completed  solution  is  fathomed  by  virtue  of  becoming 
the  new  incumbent.  Thus,  the  scheme  of  Section  4  is  finite.  I 

The  choice  of  a  search  constraint  (t)  is  the  subject  of  Section  7.  However,  before 
searching  any  of  the  constraints  binding  at  x  for  their  1-ceiiing  points,  a  preliminary 
search  is  made  for  the  1-ceiling  points  with  respect  to  a  specially-constructed  constraint 
called  an  intersection  cut,  which  is  described  in  the  next  section. 
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6.  Searching  the  Intersection  Cut 


The  intersection  cuts  employed  in  the  Exact  Ceiling  Point  Algorithm  were  developed 
by  Balas  [4]  as  a  class  of  cutting  planes  for  solving  integer  programming  problems.  They 
possess  the  usual  feature  of  chopping  off  the  optimal  ( LPr  )  solution  without  cutting  off  any 
feasible  integer  solutions  for  ( ILP ).  Furthermore,  they  make  use  of  the  same  structural 
information  from  the  ( LPr )  as  do  the  Ceiling  Point  Algorithms  described  here.  Our 
interest  in  these  cuts  stems  not  from  their  ability  to  solve  (ILP)'s  as  part  of  an  iterative 
cutting  plane  procedure,  but  from  their  ability  to  define  a  region  near  x  which  is  likely  to 
contain  a  near-optimal  or  even  optimal  feasible  integer  solution.  Thus,  an  intersection  cut 
is  introduced  into  (ILP)  after  the  Heuristic  Ceiling  Point  Algorithm  has  been  executed. 
A  suitable  region  associated  with  this  intersection  cut  is  searched  for  its  1-ceiling  points 
as  a  heuristic  step  that  uses  the  framework  of  the  exact  enumeration  scheme  described  in 
Section  4. 

The  spherical  intersection  cut  introduced  in  [4]  is  derived  from  a  convex  set  denoted 
HS(UHC[x]),  the  hypersphere  circumscribing  UHC[x],  and  is  illustrated  in  Figure  6.  A 
companion  article  [5]  describes  another  intersection  cut  derived  from  a  larger  region  called 
the  dual  to  the  unit  hypercube,  denoted  DU HC[x ].  One  cut  or  the  other  is  employed  in  the 
Exact  Ceiling  Point  Algorithm,  depending  upon  the  size  of  the  problem.  The  intersection 
cut,  denoted  by  (/),  plays  a  role  similar  to  that  of  the  functional  search  constraints  de¬ 
scribed  previously.  We  first  define  a  search  hyperrectangle  for  the  intersection  cut,  denoted 
SHR(I),  and  proceed  to  enumerate  its  feasible  1-ceiling  points.  The  only  difference  from 
the  case  of  ordinary  search  constraints  is  that  our  SHR(I)  is  not  so  large  as  to  contain  all 
of  the  1-ceiling  points  with  respect  to  (/),  but  only  the  ones  most  likely  to  be  feasible. 

The  search  region  SHR(I)  depends  upon  a  set  of  points  {r*,  k  =  1,2,...}  located 
along  the  extreme  rays  of  the  feasible  region.  Each  rk  represents  the  intersection  of  the 
kth  extreme  ray  and  the  HS(U  HC[x]).  The  intersection  cut  is  then  constructed  so  as 
to  pass  through  this  set  of  points.  These  points  also  provide  an  upper  bound  z  on  the 
optimal  objective  function  value  for  (ILP):  z  =  max*{crr*}.  The  point  rk  hats  the  form 
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Figure  6.  An  Intersection  Cut  derived  from  HS(U  HC[x]). 


x  +  A kdk,  where  A*  is  defined  such  that  x  +  A kdk  lies  on  the  surface  of  H S{UHC\x\). 
Precise  expressions  for  the  A*  are  given  in  [4]  and  [5].  Each  of  these  non-integer  points  rk 
is  rounded  to  an  integer  solution  f*  (a  vertex  of  UHC[rk])  which  is  feasible  with  respect 
to  the  intersection  cut. 

The  search  hyperrectangle  SHR(I)  is  then  determined  as  a  set  of  lower  and  upper 
bounds  =  1, . . . ,n}  by  finding  the  minimum  and  maximum,  respectively,  over 

the  jtk  component  of  all  the  f*’s.  As  with  the  bounds  defining  SHR(i ),  the  bounds 
defining  SHR(I)  should  be  no  wider  than  those  defining  SHR(S)  since  searching  beyond 
the  boundary  of  SHR(S)  cannot  improve  upon  the  incumbent.  More  precisely, 

=  min{min{f  *},£;},  iocj  = 

and 

7 >  for  j  =  1, . . .  ,n. 
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If,  for  any  component  j,  we  find  7j  <  I_j,  then  SHR{I)  =  0,  *.e.,  there  are  no  integer 
solutions  at  all  in  this  search  hyperrectangle.  In  this  case,  the  exact  algorithm  proceeds  to 
search  for  1-ceiling  points  with  respect  to  an  original  functional  constraint.  Assuming  this  is 
not  the  case,  we  enumerate  the  integer  solutions  contained  within  the  search  hyperrectangle 
SHR(I )  in  the  same  manner  as  described  in  Section  4  for  an  ordinary  search  constraint 

(*)• 

A  nice  feature  about  the  intersection  cut  is  that  the  volume  of  SHR(I)  reflects  the 
shape  of  the  feasible  region  near  x.  When  the  cone  Fit  with  vertex  x  is  narrow,  the  search 
hyperrectangle  SHR(I)  is  apt  to  contain  a  relatively  small  number  of  integer  solutions 
(feasible  or  not)  and  so  not  much  time  will  be  wasted  searching  an  unpromising  region. 
On  the  other  hand,  when  the  cone  FR  is  fairly  wide,  SHR(I)  is  apt  to  contain  a  relatively 
large  number  of  integer  solutions  and  the  subsequent  search  will  be  thorough  in  a  part  of 
the  feasible  region  with  a  good  chance  of  containing  an  optimal  solution. 


7.  Choice  of  Search  Constraint 


We  now  describe  our  rule  for  how  a  search  constraint  is  selected  in  the  Exact  Al¬ 
gorithm.  As  background,  note  that  once  constraint  (t)  has  been  chosen  and  the  “band” 
or  volume  of  the  feasible  region  between  (t)  and  (i1)  has  been  exhaustively  searched  for 
l-CP(t)’s,  constraint  (t)  may  be  replaced  by 

ofx  <bi-  max  |a,y |  ( i ") 

i 

for  purposes  of  continuing  the  search  elsewhere  in  the  feasible  region.  Therefore,  when 
this  stage  in  the  algorithm  is  reached,  ( LPr )  should  be  resolved  with  6,  replaced  by  b,"  — 
bi  —  maxy  |aj_j|,  leading  to  a  reduction  in  the  optimal  objective  function  value  of  (LPr). 
Let  (LP^)i  denote  this  new  but  related  (LPr)  whose  optimal  solution  value  is  z\.  Recall 
that  z  is  the  optimal  objective  function  value  of  (LPr)  and  serves  as  an  upper  bound  on 
the  optimal  objective  function  value  of  ( ILP ).  To  estimate  the  extent  to  which  z  and  z\ 
aiffer  without  actually  solving  (LP'R)i,  let  tt,  be  the  shadow  price  of  binding  constraint 
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(i)  for  (LPr),  where  this  shadow  price  measures  the  rate  at  which  the  objective  function 
z  changes  as  6j  is  changed  by  a  small  amount.  If  A =  b,  —  b,»  =  max_,  |at;|,  then 
z  —  z'i  >  itiAbi.  The  inequality  is  necessary  because  constraint  (t")  may  not  be  binding  at 
the  optimal  solution  to  ( LP'R)i .  However,  we  approximate  z\  by  z "  =  z  —  7r,  A5,. 

Now  let  A  =  {i|  af x  =  b{)  be  the  set  of  constraints  binding  at  x  and  let  A'  =  {i  E 
A\  z\  <  2},  where  2  is  the  objective  function  value  of  the  incumbent.  Thus,  A'  includes  any 
given  binding  constraint  only  if  fathoming  that  constraint  would  indicate  that  the  problem 
has  been  solved  because  2  is  at  least  as  large  as  the  objective  function  value  of  any  solution 
in  the  remaining  ( ILP ).  Our  rule  for  choosing  a  search  constraint  is  to  select  the  binding 
constraint  ( h )  such  that 

h  E  arg  min  { z"  } . 

The  motivation  for  this  rule  is  that  fathoming  constraint  ( h ),  ».e.,  exhaustively  searching 
SHR{h)  for  l-CP(h)'s,  gives  the  largest  guaranteed  reduction  in  2.  If,  in  addition,  h  E  .4', 
then  fathoming  constraint  (h)  may  entail  more  searching  than  is  absolutely  necessary  to 
solve  the  problem  as  a  result  of  the  ceiling  point  constraint  ( h ")  associated  with  (h)  cutting 
too  deeply  into  the  feasible  region.  To  prevent  any  unnecessary  searching,  the  right  hand 
side  bk"  for  ceiling  point  constraint  ( h ")  is  modified  (assuming  that  all  Cj  are  integer)  so 
that  z'l  =  2  +  1: 

h"  =  bk  -(*-(*  + !))/**• 

In  this  case,  we  will  be  seeking  d-CP(h)’s,  where  d  >  1.  On  any  subsequent  iterations, 
constraint  ( h ")  would  not  be  considered  for  selection  as  the  search  constraint. 

Although  similar,  this  rule  for  choosing  a  search  constraint  and  tha.  proposed  for 
the  Heuristic  Ceiling  Point  Algorithm  in  (24,  Section  3.1]  may  select  different  search  con¬ 
straints.  The  rule  presented  earlier  places  more  emphasis  on  how  the  objective  function 
value  changes  over  just  the  feasible  portion  of  the  constraint.  Here,  the  emphasis  is  on 
the  effect  of  fathoming  a  particular  constraint  and  its  associated  search  hyperrectangle 
SHR(i );  the  selected  constraint  is  that  which  guarantees  the  largest  decrease  in  2,  and 
hence,  the  largest  decrease  in  the  value  of  the  upper  bound  on  the  optimal  objective  func¬ 
tion  value  of  the  remaining  (ILP).  This  emphasis  increases  the  likelihood  that  only  a 
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small  number  of  constraints  must  be  searched  for  their  1 -ceiling  points  before  the  problem 
is  solved.  It  is  also  clear  that  the  better  the  quality  of  the  solution  that  the  Exact  Ceiling 
Point  Algorithm  starts  with,  the  fewer  iterations  it  is  likely  to  require  to  solve  the  problem. 


8.  Other  Computational  Issues 


Dining  the  computational  testing  of  the  Exact  Ceiling  Point  Algorithm  it  was  often 
noticed  that  the  same  number  of  integer  solutions  were  contained  within  the  search  hy¬ 
perrectangles  SHR(i)  and  SHR(S),  i.e.,  JJy(l  +  u}  —  lj)  =  +  ly}  —  S;),  probably 

indicating  that  the  two  hyperrectangles  coincided.  This  seemed  to  imply  that  either  the 
current  solution  was  very  close  to  being  optimal,  leading  to  a  relatively  small  SHR(S),  or 
that  the  distance  between  the  constraint  hyperplanes  (i)  and  (i1)  was  fairly  wide,  thereby 
generating  a  relatively  large  SHR(i),  or  a  combination  of  both.  Under  these  circumstances 
it  seemed  wise  to  simply  search  the  hyperrectangle  associated  with  the  n-simplex,  SHR(S), 
for  any  feasible  integer  solution  and  not  search  only  for  feasible  l-CP(i)’s.  Any  such  so¬ 
lution  found  at  this  stage  would  mostly  likely  be  a  feasible  1-ceiling  point  with  respect 
to  either  (i)  or  some  other  constraint  binding  at  x.  Furthermore,  the  algorithm  can  be 
terminated  once  all  solutions  within  the  n-simplex  have  been  enumerated. 


9.  Computational  Experience 


This  section  presents  our  computational  experience  with  the  Exact  Ceiling  Point  Al¬ 
gorithm,  using  the  relevant  parts  of  [8]  as  a  guide  to  reporting  our  results.  Having  already 
described  our  algorithmic  approach,  we  begin  in  the  next  subsection  by  describing  its 
computer-based  implementation.  Subsection  9.2  discusses  the  experimental  design,  in¬ 
cluding  the  objective  of  the  experiment,  the  origin  of  all  of  the  test  problems  used  and 
the  choice  of  performance  indicators.  Computational  results  with  the  Exact  Ceiling  Point 
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Algorithm  are  reported  in  subsection  9.3. 
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9.1.  Computer  Implementation 


Computational  testing  of  the  algorithms  was  performed  on  a  Digital  Equipment  Cor¬ 
poration  VaxStation  II  with  ten  megabytes  of  main  memory,  under  the  MicroVMS  oper¬ 
ating  system,  version  4.5.  All  of  the  code  was  written  in  Fortran  and  compiled  with  the 
VAX  Fortran  Compiler,  version  4.5,  using  the  default  settings  that  include  an  optimizer. 
Real  variables  were  declared  as  double  precision  variables.  Groups  of  test  problems  were 
submitted  as  a  batch  job  in  order  to  maintain  consistent  timing  results. 

A  clock-reading  routine  due  to  [20]  returning  CPU  time  in  centiseconds  was  employed 
to  establish  execution  times  of  various  parts  of  the  code.  Thus,  execution  times  reported 
for  the  Ceiling  Point  Algorithms  are  accurate  to  at  most  0.01  CPU  seconds.  However, 
it  is  felt  that  this  relatively  small  uncertainty  in  the  timing  can  be  safely  ignored  in  the 
following  analysis.  All  execution  times  are  given  in  CPU  seconds.  Those  execution  times 
that  apply  specifically  to  the  Exact  Ceiling  Point  Algorithm  include  the  time  required  to 
read  in  the  data  but  not  to  write  out  any  information;  those  reported  for  other  algorithms 
may  or  may  not  include  input /output  time. 


9.2.  Experimental  Design 


The  main  objective  of  our  computational  testing  was  to  assess  whether  or  not  the 
methods  for  enumerating  1-ceiling  points  described  in  this  report  constitute  a  practical 
approach  for  exactly  solving  general  integer  linear  programming  problems.  To  assess  the 
effectiveness  of  the  Exact  Ceiling  Point  Algorithm,  we  shall  compare  its  performance  to 
those  of  other  algorithms  on  a  common  set  of  test  problems.  It  should  be  emphasized  that 
these  computational  results  provide  only  a  general  indication  of  an  algorithm’s  performance 
rather  than  conclusive  evidence  because  not  only  are  we  examining  performance  based  on  a 
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relatively  limited  amount  of  computational  experience,  but  also  the  algorithms  have  been 
coded  by  different  authors,  run  on  computers  of  different  generations  and  sizes,  and  so 
forth. 

The  48  test  problems  taken  from  the  literature  have  been  grouped  into  two  categories: 
“realistic”  (because  these  problems  were  drawn  from  real  applications)  and  “randomly 
generated”  (because  the  parameters  of  these  problems  were  randomly  generated).  Charac¬ 
teristics  of  the  sets  of  realistic  and  random  test  problems  are  shown  in  Tables  1(a)  and  1(b), 
respectively.  The  first  two  columns  of  each  table  give  the  size  of  the  constraint  matrix  (rows 
by  columns)  and  the  name  by  which  we  shall  refer  to  each  problem.  Density  is  simply  the 
percentage  of  coefficients  of  the  constraint  matrix  which  are  nonzero.  A  negative  entry  in 
the  column  of  optimal  objective  function  values  indicates  that  the  problem  originally  was 
in  the  form  of  a  minimization  rather  than  a  maximization.  The  last  two  columns  provide 
two  measures  of  the  distance  between  the  optimal  objective  function  values  for  ( ILP )  and 
(LPr).  The  first  measure  is  the  normalized  duality  gap, 

D(x,x*)  =  (cTx  -  ctx*)/||c||2, 

where  ||c||2  is  the  Euclidean  norm  of  c.  This  quantity  measures  the  Euclidean  distance 
between  the  optimal  objective  function  hyperplanes  for  (ILP)  and  ( LPr ),  i.e.,  between 
cTx  =  z  and  cTx  —  z*.  It  is  a  reasonably  good  guide  for  indicating  the  difficulty  of  the 
problem:  the  larger  the  normalized  duality  gap,  generally  the  more  difficult  it  is  to  find 
an  optimal  integer  solution  and  prove  its  optimality.  The  second  measure  is  the  duality 
gap  in  percentage  terms,  100  x  (cTx  —  cTx*)/cTx,  which  provides  some  perspective  on 
the  importance  of  actually  finding  am  optimal  integer  solution  for  a  particular  problem 
once  a  good  feasible  solution  has  been  discovered.  An  integer  solution  found  to  be  within 
some  small  percentage  of  the  optimal  LP-relaxation  objective  function  value  may  be  “close 
enough”  for  all  practical  purposes. 

All  24  of  the  realistic  problems  appeared  in  the  study  by  Trauth  and  Woolsey  [27]. 
These  consist  of  ten  fixed-charge  problems,  {FC-1,  FC-2,...,  FC-10},  five  of  the  IBM  test 
problems,  (IBM-1,  IBM-2,...,  IBM-5},  and  nine  allocation  problems,  (AL-55,  AL-60,..., 
AL-100}.  The  set  of  allocation  problems  are  all  the  same  0-1  knapsack  problem  except 
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Table  1(a).  Realistic  Test  Problem  Characteristics. 


Optimal  Value  for 

Duality  Gap 

m  x  n 

Problem 

Density 

LPr.z 

ILP  :  z • 

Norm.,o) 

4x5 

FC-1 

8.79 

7 

1.032 

20.5 

4x5 

FC-2 

9.61 

8 

0.931 

16.7 

4x5 

FC-3 

11.81 

10 

1.046 

15.3 

4x5 

FC-4 

9.22 

8 

0.704 

13.0 

6x5 

FC-5 

: 

88.61 

76 

7.279 

14.2 

6x5 

FC-6 

118.13 

106 

7.003 

10.2 

4x5 

FC-7 

70 

88.61 

76 

7.279 

14.2 

4x5 

FC-8 

70 

118.13 

106 

7.003 

10.2 

6x6 

FC-9 

50 

12.00 

9 

1.732 

mm 

FC-10 

50 

18.71 

17 

0.698 

KOI 

7x7 

-8 

0.189 

6.7 

7x7 

IBM-2 

WSM 

-7 

0.472 

20.7 

3x4 

IBM-3 

100 

-187 

0.271 

4.0 

15  x  15 

IBM-4 

53 

-9.25 

-10 

0.194 

7.5 

15  x  15 

IBM-5 

53 

-12.88 

-15 

0.549 

16.3 

11  x  10 

mm 

50.30 

50 

0.008 

11  x  10 

KS 

54.50 

52 

0.063 

4.6 

11  x  10 

AL-65 

18 

58.67 

57 

0.042 

2.9 

11  x  10 

AL-70 

18 

62.83 

62 

0.021 

1.3 

11  x  10 

AL-75 

18 

67.00 

67 

0.000 

0.0 

11  x  10 

AL-80 

18 

70.60 

68 

0.065 

3.7 

11  x  10 

AL-85 

18 

74.20 

70 

0.105 

5.7 

11  x  10 

AL-90 

18 

77.80 

75 

0.070 

3.6 

11  x  10 

AL-100 

18 

85.00 

85 

0.000 

0.0 

(°)  Gives  the  normalized  duality  gap:  D(x,x*)  =  (cTx  —  cTx*)/ 1 |c| I2 ■ 
W  Gives  the  duality  gap  in  %  terms:  100  x  ( cTx  —  cTx*)/cTx. 
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Table  1(b).  Randomly  Generated  Test  Problem  Characteristics. 


Optimal  Value  for 

Duality  Gap 

m  x  n 

Problem 

Density 

LPr  :  z 

ILP  :  z* 

Norm/a) 

mam 

15  x  15 

___ _ 

100 

2956.1 

Mmmmu 

■ 

15  x  15 

B9 

100 

1 

i 

15  x  15 

in 

100 

6171 

Bf  9 

■ 

15  x  15 

1-6 

100 

2289.1 

2234 

0.333 

1 

15  x  15 

msm 

100 

1896.3 

1875 

■ mm 

i.i 

15  x  15 

100 

1758.8 

1725 

1.9 

15  x  15 

100 

2029.9 

1983 

wm 

2.3 

15  x  15 

■mm 

100 

2478.0 

2429 

wmm 

2.0 

15  x  15 

II-5 

100 

1574.8 

1558 

KSS9 

1.1 

15  x  15 

■KOI 

100 

1575.5 

1556 

1.2 

15  x  15 

BH 

100 

2088.0 

2056 

0.147 

1.5 

15  x  15 

■ 

100 

1592.8 

1548 

0.199 

2.8 

15  x  15 

n-9 

100 

1756.8 

1743 

0.063 

15  x  15 

11-10 

100 

1764.7 

1734 

0.137 

1.8 

30  x  15 

mim 

100 

1522.5 

1491 

BSSB 

2.1 

30  x  15 

■  9  « 

100 

1449.9 

1424 

■ 

1.8 

15  x  30 

100 

1811.6 

1785 

0.092 

1.5 

15  x  30 

11-14 

100 

2337.4 

2309 

0.089 

1.2 

6  x  21 

II-M 

100 

594 

0.181 

7.6 

15  x  15 

III-2 

■n 

110.7 

99 

■ 

10.6 

15  x  15 

III-3 

mm 

144.5 

130 

10.0 

15  x  15 

III-4 

50 

124.3 

92 

0.157 

27.6 

15  x  15 

III-5 

45 

119.5 

97 

0.097 

18.8 

15  x  15 

III-8 

49 

123.3 

113 

0.054 

8.4 

Gives  the  normalized  duality  gap:  D( x,x*)  =  (z  —  ^*)/||c||2. 
^  Gives  the  duality  gap  in  %  terms:  100  x  (z  —  z*)/z. 
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that  the  right  hand  side  increases  from  55  to  100.  It  should  be  noted  that  the  LP-relaxations 
associated  with  two  of  the  allocation  problems,  AL-75  and  AL-100,  possess  an  all-integer 
optimal  solution.  Thus,  AL-75  and  AL-100  are  solved  immediately  by  the  Heuristic  Ceiling 
Point  Algorithm  but  are  included  in  this  study  in  order  to  compare  our  results  more 
completely  with  those  reported  elsewhere.  The  fixed-charge  and  IBM  problems  were  first 
presented  in  [13]  and,  though  small,  sure  “hard”  to  solve  in  the  sense  that  the  optimal 
solutions  for  ( ILP )  and  ( LPr )  Eire  relatively  feu1  apart,  as  indicated  by  large  values  of 
the  normalized  duality  gap.  Characteristic  of  the  fixed-charge  problems  is  that  simple 
rounding  of  x  almost  never  yields  a  feasible  integer  solution. 

With  one  exception,  all  24  of  the  problems  with  randomly  generated  coefficients  have 
been  taken  from  Hillier’s  study  [16]  and  are  fully  specified  in  [17].  These  problems  are 

labeled  as  {1-1,  1-2,  1-5,  1-6},  {II  I,  II-2,...,  11-14}  and  {III-2 . HI-5,  III-8}.  Their  integer 

coefficients  were  generated  from  a  uniform  distribution  over  the  intervals  shown  in  the 
Table  II.  The  one  additional  problem  (labeled  “II-M”)  is  similar  to  a  Type  II  problem 
except  that  the  6j’s  are  smaller.  Originally  proposed  as  a  0-1  problem  in  [21],  II-M  was 
solved  as  a  general  integer  problem  in  [19],  as  it  is  here. 

Table  II.  Coefficient  Ranges  for  Randomly  Generated  Test  Problems. 


With  large  values  of  the  right  hand  sides  (&i’s)  and  with  constraint  matrices  which  are 
essentially  100  percent  dense,  the  Type  I  and  Type  II  problems  are  not  easy  to  solve.  The 
Type  I  problems  are  especially  tough  because  approximately  40  percent  of  their  constraint 
coefficients  are  negative,  while  the  other  60  percent  tire  positive.  The  Type  III  problems, 
on  the  other  hand,  are  not  particularly  challenging.  As  shown  in  Table  1(b),  the  normalized 
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duality  gap  for  a  typical  Type  I  problem  is  roughly  two  to  three  times  as  large  as  that  for 
an  average  Type  II  problem  which,  in  turn,  is  about  twice  that  for  a  Type  III  problem. 

For  algorithms  which  solve  the  ( LPr )  associated  with  ( ILP ),  an  alternative  to  simply 
reporting  CPU  time  is  to  examine  the  ratio  of  total  CPU  time  to  CPU  time  required  to 
solve  the  LP-relaxation.  This  ratio  gives  an  idea  of  how  much  work  is  required  by  the 
entire  algorithm  in  relation  to  a  relatively  efficient  and  dependable  algorithm  (the  simplex 
method)  used  in  the  first  phase  of  the  algorithm  to  solve  (LPr).  It  also  provides  a  crude 
basis  of  comparison  for  LP-based  algorithms  which  perhaps  have  been  coded  in  different 
languages  and/or  tested  on  different  types  of  computers.  With  this  measure,  various  al¬ 
gorithms’  execution  times  are  normalized  by  the  amount  of  time  to  solve  (LPr).  It  must 
be  emphasized,  however,  that  the  LP  solvers  embedded  within  the  respective  integer  pro¬ 
gramming  algorithms  may  have  been  designed  and  implemented  quite  differently,  causing 
such  comparisons  to  be  rather  rough. 

In  evaluating  exact  algorithms  for  (ILP),  the  most  important  performance  indicator 
would  seem  to  be  the  total  CPU  time  since  all  algorithms  being  compared  presumably 
find  an  optimal  solution.  It  may  also  be  meaningful  to  express  the  total  CPU  time  as  a 
fraction  of  the  time  required  to  solve  the  associated  (LPr).  In  contrast  to  some  other  areas 
of  mathematical  programming,  such  as  nonlinear  programming,  the  numerical  accuracy  of 
an  optimal  solution  found  by  the  Exact  Ceiling  Point  Algorithm  is  not  really  a  subject  of 
much  concern  since  it  works  only  with  integer  solutions. 


9.3.  Results  with  the  Exact  Ceiling  Point  Algorithm 


Tables  6-III(a)  and  (b)  indicate  the  relative  performance  of  the  various  phases  of  the 
Exact  Ceiling  Point  Algorithm  on  the  realistic  and  randomly  generated  test  problems, 
respectively.  The  columns  labeled  “ LPr ,”  “Heur.”,  “I-Cut,”  and  “Exact”  give  the  fraction 
of  the  total  CPU  time  spent  in  the  various  components  of  the  algorithm:  solving  the  linear 
programming  relaxation  associated  with  (ILP),  running  Phases  2  and  3  of  the  Heuristic 
Ceiling  Point  Algorithm,  searching  the  region  defined  by  the  Intersection  Cut  (see  Section 
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Table  6-III(a).  Execution  Times  of  the  Exact  Ceiling  Point  Algorithm 

on  Realistic  Problems. 


%  of  Total  CPU  time 

Total 

Problem 

LPr 

Heur. 

I-Cut 

Exact 

CPU  time 

Ratio<“> 

50.0 

10.7 

14.3 

■HE  3HE 

I 

69.6 

26.1 

4.3 

0.0 

1.4 

56.5 

13.0 

17.4 

0.23 

1.8 

FC-4 

76.5 

17.6 

5.9 

0.0 

0.17 

1.3 

FC-5 

22.1 

23.5 

2.9 

51.5 

0.68 

4.5 

FC-6 

25.4 

23.7 

3.4 

47.5 

0.59 

3.9 

FC-7 

21.7 

25.0 

3.3 

50.0 

0.60 

4.6 

FC-8 

27.3 

25.4 

1.8 

45.5 

0.55 

3.7 

FC-9 

51.5 

15.1 

6.1 

27.3 

0.33 

1.9 

FC-10 

25.0 

51.7 

14.7 

8^ _ 

1.16 

4.0 

mmm 

■ 

18.0 

■ 

■VjljTiVB 

mu 

1 

I  I 

is 

60.0 

28.0 

8.0 

IBM-4 

24.2 

75.8 

0.0 

I 

2.69 

■SB 

IBM-5 

1.3 

4.1 

6.4 

77.2 

Ifl 

0.0 

i mmm 

1 

10.7 

KM 

AL-65 

6.8 

■H 

AL-70 

56.7 

0.0 

0.0 

0.67 

AL-75 

100.0 

0.0 

0.0 

0.0 

0.28 

1.0 

39.7 

45.2 

8.2 

6.9 

0.73 

2.5 

AL-85 

39.5 

44.7 

7.9 

7.9 

0.76 

2.5 

AL-90 

38.9 

47.2 

6.9 

6.9 

0.72 

2.6 

AL-100 

100.0 

0.0 

0.0 

0.0 

0.29 

1.0 

Ratio  =  (Total  CPU  time)/(CPU  time  solving  LPr) 
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Table  6-III(b).  Execution  Times  of  the  Exact  Ceiling  Point  Algorithm 
on  Randomly  Generated  Problems. 


%  of  Total  CPU  time 

Total 

Problem 

LPr 

Heur. 

I-Cut 

Exact 

CPU  time 

Ratio 

0.1 

-fgJJk 

99.4 

348.49 

645.4 

1.1 

mam. 

85.9 

<0.1 

<0.4 

<0.1 

>99.2 

>2000 

1-6 

<0.1 

<0.2 

<0.1 

>99.6 

>900 

>1837 

II- 1 

46.5 

12.7 

2.13 

4.8 

II-2 

21.6 

27.9 

39.9 

4.6 

II-3 

7.4 

3.1 

78.9 

6.79 

13.6 

II-4 

17.5 

5.0 

67.2 

3.99 

9.7 

II-5 

28.5 

53.5 

11.1 

6.9 

1.44 

3.5 

■  LO 

1.5 

6.7 

91.2 

28.49 

66.3 

MSm 

1.8 

2.6 

1.2 

94.4 

23.35 

54.3 

II-8 

0.4 

2.2 

96.9 

116.26 

252.7 

II-9 

22.5 

49.3 

9.6 

18.7 

2.09 

4.5 

11-10 

9.0 

6.4 

54.0 

4.98 

11.1 

' 

■mi 

8.8 

37.6 

6.17 

8.8 

. 

i 

1 1 

98.8 

231.50 

<0.3 

>99.4 

>900 

>1154 

;  r 

i 

23.4 

4.4 

21.2 

30.48 

B9 

mum 

5.0 

8.0 

84.8 

13.13 

m-2 

■ 

18.3 

8.8 

2.8 

III-3 

13.2 

7.9 

mBm 

3.4 

III-4 

20.5 

12.4 

6.5 

WHOM 

4.9 

III-5 

40.0 

31.0 

15.0 

14.0 

1.00 

2.5 

III-8 

38.3 

34.0 

18.1 

9.6 

0.94 

2.6 

Ratio  =  (Total  CPU  time)/(CPU  time  solving  LPr ) 
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6),  and  executing  the  iterative  portion  of  the  Exact  Ceiling  Point  Algorithm  (see  Section 
5).  The  next  to  last  column  gives  the  total  CPU  time  in  seconds,  while  the  last  column 
gives  the  ratio  of  the  total  execution  time  to  the  amount  of  time  solving  ( LPr ). 

Based  on  both  total  CPU  time  and  the  ratio  of  total  CPU  time  to  time  spent  solving 
(LPr),  the  Exact  Ceiling  Point  Algorithm  appears  to  be  an  efficient  method  for  solving  all 
of  the  realistic  test  problems  except  IBM-5.  Aside  from  this  one  problem,  the  only  problems 
where  more  than  half  the  total  CPU  time  was  spent  outside  of  the  Heuristic  Algorithm 
were  {FC-5,...,  FC-8}.  These  four  problems  are  poorly  scaled  versions  of  {FC-1,...,  FC-4} 
in  the  sense  that  the  right  hand  sides  of  some  of  their  constraints  have  been  multiplied  by 
a  factor  of  ten,  greatly  enlarging  the  feasible  region.  Furthermore,  the  normalized  duality 
gap  for  each  of  these  problems  is  particularly  large.  In  solving  each  of  these  four  problems, 
all  four  functional  constraints  binding  at  x  had  to  be  searched  for  1-ceiling  points.  Thus, 
more  time  was  spent  in  the  iterative  portion  of  the  exact  algorithm  on  these  problems 
than  in  the  other  problems.  As  can  be  seen  in  Table  6-IV(a),  the  Intersection  Cut  search 
proved  ineffective  on  all  problems  except  FC-10,  another  indication  of  the  difficulty  of  the 
fixed-charge  test  problems. 

Table  6-III(b)  indicates  that  the  performance  of  the  Exact  Ceiling  Point  Algorithm 
on  the  randomly  generated  problems  varied  noticeably  with  the  type  of  problem.  This 
is  partly  due  to  the  fact  that  the  same  is  true  of  the  Heuristic  Ceiling  Point  Algorithm, 
which  provides  the  Exact  Algorithm  with  an  initial  feasible  integer  solution.  The  better  the 
objective  function  value  of  this  initial  solution,  the  smaller  the  size  of  the  region  examined 
by  the  Exact  Ceiling  Point  Algorithm,  and  hence,  the  quicker  it  solves  the  problem. 

On  all  of  the  Type  III  problems,  the  Heuristic  Algorithm  identified  an  optimal  solution 
and  the  Exart  Algorithm  was  able  to  prove  optimality  rather  quickly.  Of  the  fifteen  Type 
II  problems,  the  Heuristic  Algorithm  located  an  optimal  solution  in  six  cases.  On  the 
remaining  nine  Type  II  problems,  the  Intersection  Cut  search  was  effective  in  locating  a 
better  solution  four  times,  one  of  which  was  optimal.  While  only  moderately  successful, 
the  Intersection  Cut  search  is  cheap  enough  computationally  to  make  it  worthwhile.  The 
iterative  portion  of  the  Exact  Ceiling  Point  Algorithm  wets  fairly  successful  in  solving  the 
majority  of  Type  II  problems  reasonably  quickly.  However,  it  was  quite  slow  in  finishing 
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Table  6-IV(a).  Raw  Data  for  the  Exact  Ceiling  Point  Algorithm’s  Performance 

on  Realistic  Problems. 
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off  both  II-8  and  11-12  and  failed  to  prove  optimality  on  the  30-variable  problem  11-13 
within  a  900  second  execution  time  limit. 

It  should  be  noted  that  for  each  of  the  randomly  generated  problems,  a  run-time 
option  was  activated  which  reordered  the  variables  based  on  the  relative  magnitudes  of 
their  associated  cost  components.  Specifically,  the  variable  with  the  largest  Cj  was  relabeled 
as  Xi  (the  first  variable  to  be  fixed),  the  variable  with  the  second  largest  cost  coefficient 
relabeled  as  x2  (the  second  variable  to  be  fixed),  and  so  forth.  For  some  of  the  randomly 
generated  problems,  this  simple  scheme  was  very  effective  in  reducing  the  execution  time 
of  the  enumeration  scheme.  The  execution  time  of  11-12,  for  example,  was  reduced  to  one 
third  of  that  required  when  the  variables  were  not  reordered. 

The  last  set  of  tables  give  the  performance  of  the  Exact  Ceiling  Point  Algorithm 
(XCPA)  on  the  test  problems  along  with  the  performances  of  a  few  other  exact  algorithms, 
one  of  which  is  contained  within  a  widely  available  package  called  the  Generalized  Algebraic 
Modeling  System  (GAMS,  Version  2.04)  developed  by  Brooke,  Kendrick  and  Meeraus  [6]. 
When  faced  with  a  mixed  integer  linear  programming  problem,  GAMS  calls  upon  the 
Zero/One  Optimization  Methods  (ZOOM/XMP,  Version  2.0)  developed  by  Roy  Marsten. 
In  brief,  ZOOM  converts  every  (bounded)  general  integer  variable  into  a  sum  of  binary 
variables  and  applies  the  Pivot  &  Complement  heuristic  device  of  Balas  and  Martin  [7] 
to  find  an  initial  solution.  It  then  proceeds  with  an  LP-based  branch-and-bound  scheme. 
Fairly  tight  upper  bounds  on  the  variables  were  specified  in  order  to  keep  the  number  of 
binary  variables  relatively  small.  These  are  given  in  Appendix  A,  along  with  values  of  the 
GAMS/ZOOM  run-time  options  that  we  specified. 

A  blank  entry  in  a  table  indicates  that  no  time  was  reported  for  that  algorithm  on 
that  problem.  Only  the  Exact  Ceiling  Point  Algorithm  and  GAMS/ZOOM  were  executed 
on  the  same  computer  (a  VaxStation  II),  so  it  is  difficult  to  directly  compare  all  of  the 
stated  execution  times.  However,  in  one  study  [9]  of  various  computers’  performances 
in  solving  a  dense  system  of  linear  equations  using  a  standard  double  precision  Fortran 
package  (LINPACK),  the  IBM-370/168  appeared  to  be  at  least  five  times  as  fast  as  the 
IBM-370/ 158,  and  at  least  seven  times  as  fast  as  the  VAX  11/780.  Though  the  IBM- 
360/67  and  the  VaxStation  II  are  not  included  in  this  study,  a  knowledgeable  computer 
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Table  6~V(a).  Comparison  of  Performances  by  Exact  Algorithms  on  Realistic  Problems. 


XCPA 

GAMS/ZOOM 

[10] 

[2] 

[14] 

VaxStation  II 

VaxStation  II 

IBM  360/67 

370/168 

7090 

Problem 

Time 

Ratio 

Time 

Ratio 

Time 

Ratio 

Time 

Time 

FC-1 

0.28 

2.0 

4.30 

14.8 

0.32 

10.7 

mggm 

1.83 

FC-2 

0.23 

1.4 

3.13 

11.6 

0.27 

9.0 

m 

1.35 

FC-3 

0.23 

1.8 

3.41 

11.4 

0.33 

11.0 

1.88 

FC-4 

0.17 

1.3 

2.24 

7.7 

0.28 

9.3 

1.48 

FC-5 

0.68 

4.5 

8.62 

27.8 

41.77 

464.4 

MEM 

9.01 

FC-6 

0.59 

3.9 

12.68 

39.6 

22.07 

220.7 

mam 

7.57 

FC-7 

0.60 

4.6 

8.22 

24.9 

41.77 

1392.3 

7.83 

FC-8 

0.55 

3.7 

12.05 

38.9 

21.81 

727.0 

6.42 

FC-9 

0.33 

1.9 

1.92 

5.6 

0.43 

7.2 

HI 

3.23 

FC-10 

1.16 

4.0 

13.67 

14.9 

— 

9.15 

0.39 

1.9 

3.02 

6.9 

0.53 

6.6 

1.87 

0.40 

1.8 

6.46 

15.4 

0.55 

6.9 

3.02 

0.25 

1.7 

3.19 

11.8 

0.21 

10.5 

2.87 

IBM-4 

2.69 

4.1 

33.97 

22.8 

11.67 

IBM-5 

47.10 

77.2 

>900 

>671 

9.28 

15.5 

66.48 

AL-55 

■H 

3.6 

1.04 

2.7 

0.18 

2.42 

AL-60 

H  1 

1.9 

1.04 

2.7 

0.19 

5.08 

AL-65 

MM 

2.4 

1.26 

3.0 

0.14 

3.90 

AL-70 

0.67 

2.3 

1.27 

3.2 

0.18 

2.60 

AL-75 

0.28 

1.0 

0.34 

1.0 

0.15 

1.87 

AL-80 

0.73 

2.5 

1.02 

2.8 

0.15 

3.87 

AL-85 

0.76 

2.5 

1.33 

3.2 

0.17 

7.88 

AL-90 

0.72 

2.6 

1.22 

3.9 

0.16 

4.52 

AL-100 

0.29 

1.0 

0.31 

1.0 

0.17 

1.87 

Ratio  =  (Total  CPU  time)/(CPU  time  solving  LP/j) 
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Table  6-V(b).  Comparison  of  Exact  Algorithms  on  Randomly  Generated  Problems. 


XCPA 

GAMS/ZOOM 

[16] 

[2] 

VaxStation  II 

VaxStation  II 

IBM-360/67 

370/168 

Problem 

Time 

Ratio 

Time 

Ratio 

Time 

Ratio 

Time 

mm 

348.49 

645.4 

429.96 

462.3 

17.32 

19.0 

3.18 

H|| 

40.37 

85.9 

407.84 

261.4 

2.94 

3.2 

1.10 

>900 

>2000 

>900 

>612 

19.49 

20.5 

1-6 

>900 

>1837 

>900 

>967 

9.34 

7.5 

mm 

2.13 

4.8 

83.48 

132.5 

2.68 

3.9 

0.92 

19 

2.08 

4.6 

64.92 

95.5 

2.33 

2.9 

0.98 

6.79 

13.6 

79.34 

90.2 

2.21 

3.8 

1.20 

1 1-4 

3.99 

9.7 

83.32 

119.0 

2.28 

3.2 

1.02 

II-5 

1.44 

3.5 

59.30 

76.0 

1.48 

2.4 

0.81 

11-6 

28.49 

66.3 

78.39 

137.5 

3.30 

5.4 

■■ 

11-7 

23.35 

54.3 

74.97 

88.2 

3.22 

5.4 

BEfl 

11-8 

116.26 

252.7 

93.45 

118.3 

4.22 

6.2 

Mem 

11-9 

2.09 

4.5 

47.83 

83.9 

1.85 

2.9 

0.68 

11-10 

4.98 

11.1 

97.58 

112.2 

2.53 

3.1 

0.85 

msm 

6.17 

8.8 

6.48 

2.0 

lay 

231.50 

308.7 

7.17 

2.1 

■SB 

>900 

>1154 

162.25 

144.9 

11-14 

30.48 

37.2 

7.38 

6.7 

II-M 

13.13 

45.3 

mm 

HEBi 

2.8 

4.40 

6.7 

■SB 

2.4 

19 

BE 

3.4 

2.65 

4.1 

as 

2.5 

— 

II1-4 

as 

4.9 

3.80 

7.8 

1.60 

2.7 

BUI 

III-5 

1.00 

2.5 

6.50 

11.4 

1.58 

2.8 

El 

III-8 

0.94 

2.6 

5.92 

9.6 

1.41 

2.3 

■■ 

Ratio  =  (Total  CPU  time)/(CPU  time  solving  LP/j) 
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scientist  informed  us  that  the  IBM-370/158  is  at  least  as  fast  as  the  IBM-360/67,  and  that 
the  VAX  11/780  is  directly  comparable  to  the  VaxStation  II  [25]. 

For  the  realistic  test  problems,  Table  6-V(a)  compares  the  Exact  Ceiling  Point  Al¬ 
gorithm,  the  GAMS/ZOOM  package,  Faaland  and  Hillier’s  Accelerated  Bound-and-Scan 
Algorithm  [10]  run  on  an  IBM-360/67,  Austin  and  Ghandforoush’s  Surrogate  Cutting 
Plane  Algorithm  [2]  (an  extension  of  their  Reduced  Advanced  Start  Algorithm  [1])  run 
on  an  IBM-370/168  and  Haldi  and  Isaacson’s  cutting  plane  code  LIP-1  [14]  run  on  an 
IBM- 7090.  (Reasonably  good  times  on  a  subset  of  these  problems  are  also  reported  by 
[12],  but  since  they  do  not  report  anything  for  the  more  difficult  problems,  their  limited 
results  were  not  included  here.)  Based  on  the  execution  times  and  relative  speeds  of  the 
computers,  the  Exact  Ceiling  Point  Algorithm  appears  to  be  highly  competitive  with  each 
of  the  other  algorithms  on  all  of  the  realistic  problems  except  IBM-5.  Run  on  the  fastest 
computer  of  any  shown  here,  the  Surrogate  Cutting  Plane  Algorithm’s  execution  times  are 
the  smallest  and  most  stable.  It  is  interesting  that  the  Surrogate  Cutting  Plane  Algorithm 
did  not  find  the  poorly-scaled  {FC-5,...,  FC-8}  problems  any  more  difficult  to  solve  than 
{FC-1,...,  FC-4}.  All  of  the  other  algorithms  found  the  set  of  problems  {FC-5,...,  FC-8, 
FC-10,  IBM-4,  IBM-5)  relatively  difficult  to  solve. 

For  the  randomly  generated  test  problems,  Table  6-V(b)  compares  the  Exact  Ceil¬ 
ing  Point  Algorithm,  GAMS/ZOOM,  Hillier’s  Bound-and-Scan  Algorithm  [16]  run  on  an 
IBM-360/67  and  the  Surrogate  Cutting  Plane  Algorithm.  The  total  time  figures  for  the 
Bound-and-Scan  Algorithm  were  compiled  by  adding  the  times  for  the  heuristic  algorithm 
described  in  [15]  to  those  times  given  in  [16]  since  the  best  solution  yielded  by  the  former 
was  used  as  an  initial  solution  for  the  latter  algorithm  [16,  p.  668].  The  execution  times 
reported  for  the  algorithms  of  [16]  and  [2]  are  very  stable.  However,  since  the  IBM-360/67 
is  probably  several  times  slower  than  the  IBM-370/168,  Hillier’s  algorithm  may  be  faster 
than  the  algorithm  of  Austin  and  Ghandforoush.  Overall,  both  of  these  algorithms  are 
more  consistent  than  the  Exact  Ceiling  Point  Algorithm  (and  GAMS/ZOOM)  on  problems 
of  Type  I  and  II.  While  the  Exact  Ceiling  Point  Algorithm  is  competitive  with  these  two 
algorithms  on  about  half  of  the  Type  II  problems,  GAMS/ZOOM  is  not  really  competitive 
at  all,  despite  the  fact  that  it  was  only  required  to  find  a  solution  within  2  percent  of  op- 
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timality  for  Type  II  problems.  On  Type  III  problems,  the  Exact  Ceiling  Point  Algorithm, 
the  Bound-and-Scan  Algorithm  and  the  Surrogate  Cutting  Plane  Algorithm  all  appear  to 
be  of  roughly  comparable  efficiency.  Even  on  the  0-1  problems  (classes  AL  and  Type  III) 
for  which  it  was  designed,  GAMS/ZOOM  did  not  perform  as  well  as  the  Exact  Ceiling 
Point  Algorithm,  which  was  not  designed  for  0-1  problems. 

Tables  6-V(a)  and  (b)  are  summarized  by  problem  class  in  Table  6-V(c). 

Table  6-V(c).  Summary  of  Performances  by  Exact  Algorithms. 


XCPA 

GAMS/ZOOM 

[10] 

[2] 

[14] 

VaxStation  II 

VaxStation  II 

IBM  360/67 

370/168 

7090 

Class 

Time 

Ratio 

Time 

Ratio 

Time 

Ratio 

Time 

Time 

FC 

■SI 

2.9 

7.02 

9.7 

mi 

285.2 

■a 

4.98 

IBM 

Bfifl 

17.3 

>189 

>145 

9.9 

I 

17.18 

AL 

wBm 

2.2 

1.17 

3.1 

MB 

Bn 

4.29 

XCPA 

GAMS/ZOOM 

[16] 

[2] 

VaxStation  II 

VaxStation  II 

IBM  360/67 

370/168 

Class 

Time 

Ratio 

Time 

Ratio 

Time 

Ratio 

Time 

I 

>546 

>1142 

>659 

>575 

B11PS1P 

12.6 

2.14 

II* 

19.16 

42.5 

76.26 

105.3 

3.9 

1.15 

III 

1.34 

3.2 

4.67 

7.9 

1.50 

2.5 

0.46 

*  Averages  over 


10.  Summary 


In  this  report,  we  have  presented  an  exact  algorithm  for  solving  ( ILP ):  given  enough 
computer  time,  it  could  in  theory  solve  to  optimality  any  (ILP)  for  which  x  exists  and 
is  unique.  The  Exact  Ceiling  Point  Algorithm  is  based  upon  theorems  implying  that  to 
solve  (ILP)  it  is  sufficient  to  enumerate  only  the  feasible  1-ceiling  points  lying  within  the 
n-simplex  S.  This  is  done  by  systematically  searching  for  feasible  1-ceiling  points  with 
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respect  to  one  constraint  at  a  time.  While  enumerating  such  solutions,  three  features 
described  in  Section  4  help  to  accelerate  the  fathoming  process,  including  the  construction 
of  conditional  variable  bounds  and  an  ability  to  double  backtrack.  Because  the  iterative 
part  of  the  Exact  Ceiling  Point  Algorithm  benefits  greatly  by  having  as  good  an  initial 
solution  as  possible,  a  heuristic  device  is  employed  prior  to  the  execution  of  the  iterative 
procedure.  This  device  is  to  search  for  1-ceiling  points  with  respect  to  an  intersection  cut 
located  close  to  x,  thereby  examining  the  most  promising  part  of  the  feasible  region  first. 

Our  computational  experience  leads  us  to  believe  that  the  Exact  Ceiling  Point  Al¬ 
gorithm  is  a  good  approach  for  tackling  some  problems,  but  not  the  best  approach  for 
tackling  others.  Specifically,  the  Exact  Ceiling  Point  Algorithm  performed  very  well  on 
the  set  of  small  realistic  problems,  but  did  not  perform  as  well  as  several  other  algorithms 
on  some  of  the  more  difficult  randomly  generated  problems.  On  the  set  of  realistic  prob¬ 
lems,  particularly  the  fixed-charge  problems,  the  double  backtracking  feature  (see  Section 
4.2)  of  the  exact  enumeration  scheme  seems  to  be  very  helpful  in  quickly  fathoming  large 
numbers  of  solutions.  The  variables  also  appear  to  be  arranged  in  an  order  that  is  advanta¬ 
geous  to  this  scheme.  However,  on  the  randomly  generated  problems,  the  ability  to  double 
backtrack  appears  to  be  of  little  help  because  there  is  no  pattern  of  ratios  from 

one  column  to  the  next.  Also,  when  the  Heuristic  Ceiling  Point  Algorithm  did  not  find 
a  very  good  solution,  as  was  the  case  for  most  of  the  Type  I  problems,  it  often  led  to  an 
inefficient  performance  by  the  Exact  Ceiling  Point  Algorithm.  Further  research  is  needed 
to  better  determine  what  types  of  problems  are  best  and  least  suited  to  our  ceiling  point 
approach. 
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In  order  for  GAMS/ZOOM  to  convert  each  general  integer  variable  into  a  sum 
of  binary  variables,  a  reasonably  tight  upper  bound  was  specified  for  each  general 
integer  variable,  as  shown  in  Table  VI.  The  number  n'  of  binary  variables  in  the 
transformed  problem  is  given  in  the  last  column. 


Table  VI.  Upper  Bounds  Specified  in  the  GAMS/ZOOM  Runs. 


Problem/Class 

Upper  Bounds 

n' 

FC-l,...,FC-4 

Xj  <  1,  j  =  1, 2;  Xj  <  10,  j  =  3, ...,  5 

14 

FC-5,...,FC-8 

Xj  <  1,  j  =  1,2;  Xj  <  100,  j  =  3,...,  5 

23 

FC-9 

Xj  <  1,  j  =  Xj  <  10,  ji=4,...,6 

15 

FC-10 

Xj  <  1,  j  =  1, ...» 6;  Xj  <  15,  j  =  7,...,  12 

30 

IBM-1,  IBM-2 

Xj  <7,  j  =  1, ...,  7 

21 

IBM-3 

Xj  <  31,  j  =  1, ...,  4 

IBM-4,  IBM-5 

Xj  <  3,  j  —  1, ...,  15 

30 

AL 

Xj  <  1,  j  =  1 . 10 

10 

Types  I  &  II* 

Xj  <  31,  j  =  1, ...,  15 

75 

1-5 

Xj  <  50,  j  =  1, ...,  15 

90 

Type  III 

Xj  <  1,  j  =  1, ...,  15 

15 

*  except  1-5 


GAMS/ZOOM  Options  specified  in  every  Program  file  “PROBLEM. GMS” 
OPTCA  =  0.0 

OPTCR  =  0.001  (0.020  for  all  Type  II  problems) 

GAMS/ZOOM  Options  listed  in  Specs  file  “GAMSZOOM.SPC” 

BRANCH  =  YES 
DIVE  =  YES 
EXPAND  =  3 
HEURISTIC  =  YES 

INCUMBENT  =  -1000  (+1000  for  IBM  problems) 

MAX  SAVE  =  5 
PRINT  CONTINUOUS  =  0 
PRINT  HEURISTIC  =  0 
PRINT  BRANCH  =  0  ^ 

PRINT  TOUR  =  0 
QUIT  =  NO 
SELECT  =  2 

All  other  options  assumed  their  default  values.  A  preliminary  run  was  made  for 
each  test  problem  with  PRINT  HEURISTIC  =  1  in  order  to  find  zh  =  cTifj. 


Appendix  B:  Listing  of  Fortran  Code  Implementation 


Listing  of  Fortran  Code  Implementation  of 
the  Exact  Ceiling  Point  Algorithm 
for  General  Integer  Linear  Programming 


C  PLIST.FOR:  mm(nn)  =  maximum  value  of  M(N)  Put  into  all  routines 
IMPLICIT  REAL* 8  <A-H,0-Z) 

parameter  (mm“37,  nn=31,  tol=2 . 10734D-08,  bigm=l.D5, 

-  maxit=75,  one=l.D0,  zero=0.D0) 

C  C0MLP1.F0R:  Global  vars.  created/used  in  LP  routines 
dimension  A  (mm,  nn) ,  B  (mm) ,  C  (nn) ,  INDCT  (mm) , 

-  ABAR (mm, nn+mm+mm+1) ,ABAL(nn,nn) , 

ibasis (mm) , nonbas (nn+mm) , indbas (nn+mm+mm) , initbv (mm) , 

NPIV (2) , XBAR (nn+mm+mm) ,X0 (nn) , signac (mm,nn) , 
dirs (nn,nn) ,CTXPT (nn,mm) , BFCXO (mm) .PRICES (mm) 
common/clpl/A, B,C, INDCT, INDOBJ,M,N, 

ABAR, ABAL, ibasis, nonbas , indbas, NPIV, XBAR, XO, ZO, IR, IS, 
signac, nx, dirs, ctxpt, BFCXO, PRICES 
logical*2  indbas, signac,  ctxpt,BFCXO 

C  COMHRUN.FOR:  Global  vars.  used  in  HRUN  routines 

dimension  FIS3 (mm, 1+nn) , RATES (nn) , SDIR (nn) ,CTVAL (mm) 
common /ch run /f feas, FIS 3, pet 01, RATES, SDIR, nep, CTVAL 
logical*2  ffeas,ncp 

C  COMRUN1 . FOR :  Global  vars.  primarily  used  in  RUN 

real *8  RTIMES (-1 : 30) , RPCTS (-1 : 30) , RVALS (-1 : 30) , 

AIMIN (mm) , AIMAX (mm) 

intege  r  IXSTAR ( nn ) , CORDER ( nn ) , VARORD ( nn ) 
common/crunl/IXSTAR, ZSTAR, AIMIN, AIMAX, ZUP, ALL, 

RTIMES , RPCTS , RVALS , CORDER, VARORD 

C  COMXRUN2 .FOR:  Global  vars.  used  in  XRUN  and  RUNCUT 

Dimension  LAMDA(nn),  ALPHA(nn,nn) ,  P(nn,nn),  SUB(nn) 

Dimension  LOBD (nn) , UPBD (nn) , LONS (nn) ,UPNS (nn) ,CUT(nn+l) 
common/cxrun2/LAMDA, ALPHA, P, SUB, LOBD, UPBD, LONS, UPNS, CUT 
integer  LOBD, UPBD, LONS, UPNS, SUB 
real *8  LAMDA 

C  COMXRUN3 . FOR :  Global  vars.  used  in  XRUN 
integer  irange (nn) , icase 

real*8  CMPMIN (mm, nn+1) ,RA(mm,nn) , AXINC (mm, nn) , 

LL(mm+l,nn) ,UU (mm+l,nn) , sizlim, callim, xcalls (nn) 
common /cxrun3/ irange, icase, LL,UU, CMPMIN, RA, AXINC, 
sizlim, callim, xcalls 

C  COMXRUN4 .FOR:  More  XRUN  global  vars.,  especially  XCP-related  routines 
integer  LJ, first (nn) , final (nn) , inc  (nn) , newz (nn) , 
loc (nn) ,upc (nn) , ISN(nn) 
real *8  gap (mm, nn+1) 

common/cxrun4/gap, LJ, first, final, inc,nowz, loc,upc, 
ISN,minarg,maxarg 

C  COMP RT. FOR:  Print  Switches 

common/cprint/hprint,xprint 
integer*2  hprint (25) , xprint (25) 

C  COMIO .FOR:  I/O  files 

common/cinout/inf ile, iorun, iohrun, iocut,  ioxrun,  iolp 
character*64  infile, iorun, iohrun, iocut, ioxrun, iolp 
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c - 

c  By  Robert  M.  Saltzman 
c 

c  Applies  the  Exact  Ceiling  Point  Algorithm  to  each  problem 
c  listed  in  the  file  ILPDATA.DAT. 
c 

c  Parenthetical  comments  with  Section  numbers  refer  to  parts  of: 
c  Saltzman,  R.,  "Ceiling  Point  Algorithms  for  General  Integer 
c  Linear  Programming,"  unpublished  Ph.D.  dissertation.  Dept,  of 
c  Operations  Research,  Stanford  University,  Stanford,  Calif . , 
c  December  1988. 

- - 

*• 

c 

include  1 SDISK2 : [SALTZ . ILP1] plist . f or • 
include  1 SDISK2 : [ SALTZ. ILPlJcomio. for' 
include  '$DISK2: [SALTZ. ILP1] comprt . for* 
include  '$DISK2: [SALTZ.ILP13comxrun3.for1 

open (  2, file='$DISK2: [SALTZ.ILP13ilpdata.dat',  status- • old’ ) 
open(  5, file-'$DISK2 : [SALTZ.ILPlJouthrun.dat*,  status- ’unknown’ ) 
open(  6, file-'$DlSK2: [SALTZ.ILPlJoutrun.dat',  status- ' unknown ' ) 
open (22, file- ' SDISK2 : [SALTZ . ILP1] switches .dat ' , status- 'old' ) 
c 

c  ftead  in  print  switches  and  run-time  options 

read(22, *)  <hprint ( j) , j-1, 25) 
read(22,*J  (xprint ( j) , j-1,25) 
read(22,*)  sizlim, callim 
write(*,*)  'sizlim, callim  ', sizlim, callim 
write (6, *) 'Exact  Ceiling  Point  Algorithm  Summary* 
write (6,*)'  • 

write (6, *) 'Problem  LP  ZO  Set+Heur  Z  ’, 

'  IntCut  Z  Exact  Z  Total  Ratio' 
write (6,  *)  ' -  —  -  -  ', 

c 

c  Headlines  for  Heuristic  Summary  report,  if  desired 

if  (hprint(19)  .eq.  1)  then 

write (5, *) 'Heuristic  Ceiling  Point  Algorithm  Summary' 
write (5, *)  '  ' 

write (5,*) 'Problem  LP  ZO  Set+Ph.l  Z  ', 

'  Phase2  Z  Phase3  Z  Total  Ratio' 
write (5,  *)  ' -  —  — -  -  ', 

endif 

c 

c  Loop  through  all  problems  specified  in  ILPDATA.DAT 

do  8100  ip  -  1,  40 

c  Get  the  name  of  the  input  data  file,  e.g.,  HALDI5.DAT 

read (2, 8105)  infile 
if  (infile(l:3)  .eq.  'end')  goto  8150 
if  (inf ile (1 : 3)  .eq.  'END')  goto  8150 
write  (*,*)  '*****  Starting  problem  *****  ’, infile 
c 

call  RUN 
c 

write  (*,*)  '*****  Finished  problem  *****  ', infile 
8100  continue 

8105  format (All) 

8150  continue 

stop 
end 
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c 

c 


subroutine  RUN 


c 

c  Runs  Heuristic  and  Exact  Ceiling  Point  Algorithms  for  1  problem, 
c  (Overview  of  entire  algorithm  given  in  Section  5.5.) 
c  Note:  Code  for  Heuristic  Ceiling  Point  Algorithm  subroutines 
c  (LP SOLVE,  SETUP,  HRUN)  is  listed  in  Saltzman  and  Hillier  [24] . 

c - 

c 

include  '$DISK2: [SALT2.ILPl]plist.for’ 
include  ’ SDISK2 : (SALT2 . ILP1) comio . for ' 
include  *$DISK2: [SALT2 . ILPl)comlpl .for* 
include  'SDISK2: [SALT2.ILPl]comprt.for' 
include  '$DISK2: [SALT2.ILPl]comrunl .for1 
include  '$DISK2: [SALT2 . ILP1] comxrun4 . for  * 
integer  izz(5) 
c 

open  (4, f ile-' $DISK2 : [SALT2 . ILPlloutlp.dat ' ,  status-' unknown' ) 

open (8, f ile-' $DISK2 : [SALT2 . ILPlJoutruncut.dat ' , status- 'unknown' ) 
open (9, f ile-' $DISK2 : [SALT2 . ILP1] outxrun .dat ' ,  3tatus-' unknown' ) 

c 

c  Initialize  clock  reading  routines  and  summary  values 

call  XTIMER (0, 0,0) 
do  8001  i  -  -1,  30 
rtimes(i)  -  zero 
rpcts(i)  =  zero 
rvals(i)  -  zero 
8001  continue 

c 

c  Solve  LP-relaxation  of  (ILP) 

call  XTIMER(1, 0, zero) 
call  LPSOLVE 

call  XTIMER(-1, 0, RTIMES (-1) ) 
if  (xprint  (8)  .eq.0) 

write (*,*)'***>  LPSOLV: ’, RTIMES (-1) , *  Z0=',Z0 
c 

c  Check  for  all-integer  LP  solution 

do  8004  j  -  1,  n 

if  (DABS (X0 ( j ) -IRNDWN (X0 ( j ) ) )  .gt.  tol)  goto  8008 
8004  continue 

write (6,*)' - >  LP  solution  is  all-integer  (X*-X0) ' 

do  8006  j  -  1,  n 

ixstar(j)  -  IRNDWN (X0 ( j ) ) 

8006  continue 

goto  8090 
c 

c  Initialize  (ALL,  X*,  2UP( 

8008  call  XTIMER (1, 0, zero) 

ALL  =  one 

c  ALL  =  1  =>  only  search  for  solns.  strictly  better  than  incumbent 

ZUP  -  DBLE ( IRNDWN ( Z0 ) ) 
do  8010  j  -  1,  n 
IXSTAR(j)  -  -1 
8010  continue 

c 
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c  Calculate  global  vara.  (SIGNAC,  ABAL,  DIRS,  CTXPT,  BFCXO) 

call  SETUP 

call  XTIMER(-1, 0, rtimes (0) ) 

if  (xprint(8  ).eq.O)  write(*,*>  '***>  SETUP:  ',rtimes(0) 
c 

c  Set  default  values  for  Z *  in  case  HRUN  fails  or  is  bypassed 

if  ( indob j  .eq.  1)  ZSTAR  -  zero 
if  ( indob j  .eq.  -1)  ZSTAR  -  -2.*DABS(Z0) 
if  (xprint(21)  .eq.  1)  goto  8025 
if  (xprint(22)  .eq.  1)  goto  8035 
c 

c. . 

c  Run  Heuristic  Ceiling  Point  Algorithm 

call  XTIMER( 1,0, zero) 
c 

call  HRUN 
c 

if  (ixstar(l)  .It.  0) 

-  write (*, *) ' ***>  HRUN  failed.  Initial  Z*  -', ZSTAR 
call  XTIMER(-1, 0, rtimes (1) ) 
c 

izz (1)  -  IRNDWN ( ZSTAR) 
if  (xprint (8) .eq. 0) 

write (*,*)•***>  HRUN:  rtimes (1) , 1  Z*-*,izz(l) 

c 

c  Write  out  summary  information  for  the  Heuristic  Algorithm 

if  (hprint(19)  .eq.  1)  then 

rtimes (1)  -  rtimes (11) +rtimes (12) +rtimes (13) -rtimes (0) 
write  (5, 3015)  infile, rtimes (-1) , Z0, rtimes (11) , rtimes (12) , 
rvals (12) , rtimes (13) , rvals (13) , rtimes (14) , rpcts (15) 

8015  format (IX, All, F5 .2, F7 . 1, F7 .2, 6X,F8 . 2,F6 . 0, F7 .2,F6 . 0, 2F7 .2) 

write (5, 8020)  infile, rpcts (10) , rpcts (11) , rpcts (12) , rpcts  (13) 
8020  format (IX, All, F5 .2, 6X, F8 . 2,  6X,F8 .2, 5X,  F8 .2) 

endif 
c 

if  (DBLE (izz (1) )  .ge.  ZUP)  then 

write(6,*)' - >  Heuristic  alg.  found  optimal  solution' 

izz(2)  -  izz(l) 
izz(3)  -  izz(l) 
goto  8090 
endif 
c 

c . 

c  Search  the  Intersection  Cut  prior  to  full  enumeration 

8025  IF  (XPRINT (18)  .EQ.  1)  GOTO  8090 
call  XTIMER(1, 0, zero) 
call  RUNCUT 

call  XTIMER(-1, 0, rtimes (2) ) 

C 

izz  (2)  -  IRNDWN (ZSTAR) 
if  (xprint (8) .eq.0) 

write (*,*)•***>  RUNCUT:', rtimes (2),'  Z*-',izz(2) 
if  (DBLE (izz (2) )  .ge.  ZUP)  then 

write(6,*)  ' - >  I-Cut  search  found  optimal  solution' 

izz(3)  -  izz(2) 
goto  8090 
endif 
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c 

IF  (XPRINT (19)  .EQ.  1)  GOTO  8090 
c 

c . . . 

c  Eliminate  redundant  constrainta 

call  XTIMERd, 0, zero) 
call  REDCTS 
c 

c  Run  Exact  Ceiling  Point  Algorithm 

8035  call  XRUN 

call  XTIMER(-1, 0, rtimes (3) ) 
izz  (3)  -  IRNDWN (ZSTAR) 
if  (xprint (8) .eq. 0) 

write (*,*)■***>  XRUN:  • , rtimea (3) , •  Z*-',izz(3) 
c 

c  Write  out  an  optimal  solution  &  execution  times 

8090  call  REPRT1 (izz) 
c 

return 

end 

c - 

c 

subroutine  REPRTl(iz) 

c 

c  Called  by  RUN  to  report  a  few  important  pieces  of  information 

c - 

c 

include  '$DISK2: [SALTZ.ILPl]plist.for' 
include  '$DISK2: (SALTZ.ILPlJcomio.for’ 
include  *$DISK2: (SALTZ . ILPlJcomprt . for • 
include  ’$DISK2: [SALTZ . ILPl]comlpl . for ’ 
include  1 $DISK2: [SALTZ . ILPljcororunl . for ’ 
c 

integer  ixx (nn) , ixxx (nn) , iz (5) 
c  ixx  &  ixxx  are  reordered  versions  of  IXSTAR 

c  iz  is  a  vector  of  objective  function  values 

c 

c  Put  X*  in  original  variable  order  before  writing  out 

do  8091  j  -  1,  n 

ixx(varord( j) )  =  IXSTAR (j) 

8091  continue 

do  8092  j  =  1,  n 

ixxx(corder ( j) )  »  ixx(j) 

8092  continue 
c 

write  <9,*)  ’  • 
write  (9,*)  infile 
write  (9,*)  •  Z*  »',iz<3) 
do  8093  j  -  1,  n 

write (9,*)  '  X*(',j,')  -’,ixxx(j) 

8093  continue 
c 


49 


I 


c . . . 

c  Write  out  summary  information  for  the  Exact  Algorithm 

do  8095  i  -  -1,  3 

rtimes(4)  -  rtimes(4)  +  rtimes(i) 

8095  continue 

if  (xprint (8) .eq.O) 

write <*,*) •***>  TOTAL  TIME:  ', rtimes (4) 
c 

c  Combine  times  of  SETUP  t  HRUN 

rtimes (1)  -  rtimes(O)  +  rtimes(l) 
c 

c  Get  each  part's  percentage  of  TOTAL  time  £  TOTAL/LP  ratio 

do  8096  i  -  -1,  3 

rpcts(i)  -  100 .D0*rtimes (i) /rtimes (4) 

8096  continue 

tot21p  -  rtimes (4) /rtimes (-1) 
c 

write (6, 8097}  infile, rtimes (-1) , Z0, rtimes (1) , 

-  iz  (1) , rtimes (2) , iz (2) , rtimes (3) , iz (3) , rtimes (4) ,tot21p 

8097  format (IX, All, P5 .2,F7.1,F8.2,I5,F8.2,I5,F8.2,I5, 2F8 .2) 

C 

write (6, 8098)  infile, rpcts (-1) , rpcts (1) , rpcts (2) , rpcts (3) 

8098  format (IX, All, F5 . 2, 7X,F8 . 2, 5X, F8 .2, 5X, F8 . 2) 

c 

return 

end 
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c 

subroutine  RUNOUT 
c 

c  Called  by  RUN  to  search  intersection  cut  ("I-Cut") 

c - - - 

include  'SDISK2: [SALTZ . ILP1 Jplist . for  * 
include  '$DISK2: [SALTZ . ILP1) comprt . for ' 
include  *$DISK2: [ SALTZ . ILP1] comlpl . for ’ 
include  '$DISK2: [SALTZ . ILP1) comrunl . for ' 
c 

c  Perform  preparatory  routines 

call  PRECUT 
c 

if  (zstar  .ge.  zup)  goto  4020 
ihp  »  m-1 

if  (xprint  (1)  .eq.  1)  write(8,*)  '  I-CUT 1 ,  (A  (ihp,  j) ,  j*»l,  n) 
c 

c  Find  RA  once  (since  there's  no  variable  reordering  here) 

call  ARATIO 
c 

c  Enumerate  feasible  l-CP's  with  respect  to  1-Cut 

call  XCP(ihp) 
c 

c  AIMAX ( i )  =  largest  coefficient  in  absolute  value  of 

c  constraint  (i)  over  the  support  of  C. 

do  4015  i  **  1,  m 
rmax  =  -bigm 
do  4010  j  =  1,  n 

if  (C ( j)  .eq.  zero)  goto  4010 

if  (DABS  (a  (i,  j) )  .gt.  max)  max  *  DABS  (A(i,  j) ) 

4010  continue 

AIMAX (i)  =  rmax  -  one 
4015  continue 

c 

c  Adjust  AIMAX  for  I-cut  since  its  coefficients  are  non-integer 

AIMAX (m-1)  =  AIMAX (m-1)  +  one 

c 

4020  continue 

return 
end 

- - 

c 

subroutine  PRECUT 
c 

c  Called  by  RUNOUT  to  prepare  for  intersection  cut  search 

c - 

include  ' SDISK2 : [SALTZ . ILP1] plist . for ' 
include  ' SDISK2 : [SALTZ . ILP1] comprt .for' 
include  '$DISK2: [SALTZ . ILP 11 comlpl . for ' 
include  '$DISK2: [SALTZ . ILP1] comrunl . for ' 
include  ' $DISK2 : [SALTZ . ILP1 ] comxrun2 . for ' 
c 

c  all  =  1  =>  search  only  for  strictly  better  solutions 

all  *  one 

zlevel  -  zstar+all 
if  (xprint  (2)  .eq.  1)  write  (8,  *) 

'PRECUT:  zlevel=' , zlevel, '  zup“',zup 

c 
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c  Find  (P(k)'s)  =  Set  of  intersection  points 

call  OBJCUTfz level) 
c 

if  (n  .ge.  15)  goto  4025 
if  (xprint(20)  .eq.  1)  goto  4025 
c 

c  Either  get  spherical  intersection  cut  points  (alpha) 

call  CUTPT 
goto  4030 
c 

c  or,  get  dual-UHC  intersection  cut  points  (alpha*) 

4025  call  CUTPT2 

c 

4030  continue 

c 

c  Create  SUB  «  Simple  Upper  Bounds 

call  BOUNDS 

if  (ISRSIM(m)  .eq.  1)  goto  4035 

if  (xprint (2) .eq. 1)  write(8,*)  'Shr(ns)=0  =>  Z*  =  ',zstar 
zup  =  -bigm 
goto  4090 
c 

4035  continue 

c 

c  Find  ALPHA (k)  with  largest  objective  function  value 

c  It  becomes  a  new  (smaller)  upper  bound  on  Z* 

vmax  =  -bigm 
do  4040  k  -  1,  n 

value  =  VDOT (n, C, ALPHA (1, k) ) 
if  (value  .gt.  vmax)  vmax  -  value 
4040  continue 

ZUP  =  DBLE  < IRNDWN (vmax) ) 
if  (zstar  .It.  vmax)  goto  4045 

if  (xprint (2) .eq. 1)  write  (8,*)  'Z*  >«  zup  =>  Z*=', zstar 
c  goto  4090 

c 

4045  continue 

c  Determine  the  I-cut  coefficients 

call  CUTHP 

c 

c  Append  I-Cut  &  Obj.fn.  to  bottom  of  A  matrix  and  B  vectors 

do  4050  j  *  1,  n 

A(m+l,j)  =  CUT ( j ) 

A(m+2, j)  =  -C(j) 

4050  continue 

B (m+1)  =  CUT (n+1) 

B(m+2)  *  zero 

c . 

M  =  M  +  2 

c . 

4090  continue 

return 
end 
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c 

subroutine  OBJCUT (ob jval) 
c 

c  Called  by  PRECUT  &  INCMOD  to  find  the  set  of  P(lc)  's,  where 
c  P(k)  -  point  where  the  kth  extreme  ray  intersects  objective 
c  function  hyperplane:  cx  “  objval.  (See  Section  5.2) 

c - 

include  '$DISK2: [SALT2 . ILPl]plist . for ' 
include  *$DISK2: [SALT2 . ILP1) comprt . for ' 
include  ' $DISK2 : [SALTZ . I LP1] comlpl . for ' 
include  ' $DISK2 : [SALTZ . ILP1] comxrun2 . for ' 
c 

real*8  objval, theta 
if  (xprint (3) .eq.  1) 

-  write  (8, *)  'OBJCUT: objval  =  objval, '  ZO  -',Z0 
zgap  -  ZO  -  objval 
c 

c  For  each  extreme  direction  k,  find  intersection  point  P(k) 

do  4120  k  =  1,  nx 
theta  «  bigm 

dot  -  -VDOT(n,C,DIRS(l,k)  ) 

c  dot  =  o  =>  ext.  dir.  is  II  objective  function  hyperplane 

if  (DABS (dot)  .It.  tol)  goto  4100 
theta  =  zgap/dot 
c 

4100  do  4110  j  -  1,  n 

P(j,k)  -  X0(j)  +  theta*DIRS ( j, k) 

4110  continue 

if  (xprint (3) .eq. 1) 

write  (8,*)  'OBJCUT:  Pk=’ , (P ( j, k) , j«l,n) 

4120  continue 
return 
end 

c - 

c 

subroutine  CUTPT 
c 

c  Called  once  by  PRECUT  to  find  the  set  of  ALPHA (k) ’s,  where 
c  ALPHA (k)  =  point  where  kth  extreme  ray  intersects  I-cut  hyperplane 
c  See  paper  by  E.  Balas  [Ba71]  on  Spherical  Intersection  Cuts. 

c - 

include  '$DISK2: [SALTZ . ILP1] plist . for ' 
include  ' $DISK2 : [SALTZ . ILP1] comprt . for ' 
include  ' SDISK2 : [SALTZ . ILP1] comlpl . for ' 
include  'SDISK2:  [SALTZ  .  ILP1J comxrun2 .  fo..' ' 
c 

dimension  f (nn) , fe2  (nn) , emf (nn) 
do  4200  j  «  1,  n 

f ( j!  -  X0(j)  -  DINT (X0 ( J) ) 
fe2 ( j)  -  f ( j)  -  0.5D0 
emf ( j )  -  one  -  f ( j ) 

4200  continue 

fef  -  VDOT (n, f ,emf ) 
c 
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do  4220  k  =  1,  nx 

ak2  -  VDOT  (n,  ABAL  (1,1c),  ABAL ( 1,  k)  ) 
hk  «  VDOT (n,f e2, ABAL (l,k) ) 
lamda <k)  »  DSQRT < (hk*hk) + (fef *ak2) ) 
lamda(k)  -  (hk  +  lamda (k) ) /ak2 
do  4210  j  -  1,  n 

alpha(j,k)  -  X0(j)  -  lamda (k) *abal (j, k) 

4210  continue 

if  (xprint (4) .eq.l) write (8, *) 'Alpha' , (alpha (j, k) , j-l,n) 
4220  continue 

if  (xprint (4) .eq. 1)  write(8,*)  'Lamda ', (lamda (j) ,  j-1,  n) 

return 

end 

c - 

c 

subroutine  CUTHP 
c 

c  Called  once  by  PRECUT  to  compute  I-cut  coefficients  ("CUT") 

c - 

include  '$D1SK2: (SALTZ . ILPlJplist . for ' 
include  '$DISK2: [SALTZ . ILP1] comprt . for ’ 
include  '$DISK2: (SALTZ . ILP1] comlpl . for ' 
include  '$DISK2: (SALTZ . ILPl]comxrun2 . for ' 
c 

do  4300  k  -  1,  n 
cut(k)  «  zero 

if  (.not.  indbas(k))  cut(k)  -  -one/ lamda (k) 

4300  continue 

cut (n+1)  »  -one 
c 

do  4320  j  -  1,  n 
nb  »  nonbas(j) 
if  (nb  .le.  n)  goto  4320 
i  -  nb  -  n 

if  (i  .gt.  m)  goto  4320 
do  4310  k  -  1,  n 

cut(k)  -  cut(k)  +  (A(i,k) /lamda (j) ) 

4310  continue 

cut (n+1)  •  cut (n+1)  +  (B (i) /lamda ( j) ) 

4320  continue 

c 

c  Make  sure  that  CUT  chops  off  X0  (if  not,  mult,  by  -1) 

vlhs  =  VDOT (n, CUT, XO) 
if  (vlhs  .gt.  CUT (n+1))  goto  4340 
do  4330  j-1,  n+1 
cut(j)  -  -cut(j) 

4330  continue 

4340  if  (xprint (5) .eq.l)  write  (8,*)  ’ ICUT’ , (CUT ( j) , j-1, n+1) 

return 
end 


» 
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c - 

c 

subroutine  CUTSHR 
c 

c  Called  by  ISHR  to  calculate  SHR(I-Cut),  the  set  of  unconditional 
c  bounds  for  all  variables.  They  are  based  on  the  alpha (k, j) ' s, 
c  rounded  according  to  the  sign  of  CUT(j) . 

c - 

include  '$DISK2: [SALTZ . ILPlJplist . f or ' 
include  '$DISK2:  [SALTZ  .  ILP1]  cornprt .  for  ■ 
include  ' SDISK2 : [SALTZ . ILP1] comlpl . for ' 
include  '$DISK2: [SALTZ . ILP1] comxrun2 . for ' 
c 

c  Work  through  each  component  j 

do  4430  j  =  1,  n 

lobd(j)  -  IRNDWN (bigm) 
upbd(j)  =  0 
c 

c  Consider  all  the  ALPHA(k)’s 

do  4425  k  =  1,  n 
c 

if  (CUT ( j ) )  4405,4410,4415 
c  if  cut  (k)  <  0,  round  alpha (j,k)  up 

4405  ibd  -  IRNDUP (alpha (j, k) ) 

go  to  4420 

c  if  cut(k)  is  0,  round  alpha (j,k)  to  nearest  int 

4410  ibd  =  IDNINT (alpha (j, k) ) 

go  to  4420 

c  if  cut (k)  >  0,  round  alpha (j,k)  down 

4415  ibd  =  IRNDWN (alpha ( j, k) ) 

c 

4420  if  (ibd  .It.  lobd(j))  lobd(j)  =  ibd 

if  ( ibd  . gt .  upbd ( j ) )  upbd ( j )  =  ibd 
4425  continue 

c 

c  Take  intersection  of  n-simplex  and  uncdl.  bds . 

if  (lons(j)  .gt.  lobd(j))  lobd(j)  =  lons(j) 
if  (upns(j)  .It.  upbd(j))  lobd(j)  =  upns(j) 

4430  continue 

c 

4450  continue 

return 
end 

c - 

c 

subroutine  ARATIO 
c 

c  Called  by  PRECUT  and  REORDR  to  calculate  the  ratio  of  adjacent 
c  elements  of  A(i,),  row  by  row  (■=  f(i,j)  in  Section  5.4).  Also 
c  sets  up  global  variables  INC  and  AXINC. 

c - 

include  ' SDISK2 : [SALTZ . ILP1 ]plist . for ’ 
include  '$DISK2: [SALTZ . ILP1] cornprt . for ' 
include  ' $DISK2 : [SALTZ . ILP1] comlpl . for ’ 
include  ’$DISK2: [SALTZ . ILPl]comxrun2 . for ' 
include  'SDISK2: [SALTZ . ILP1] comxrun3 . for ' 
include  ' SDISK2 : [SALTZ . ILP1] comxrun4 .for’ 
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c 

c 

Determine  INC(j)=  direction  in  which  to  enumerate 
do  4510  j  -  1,  n 

X { j )  values 

c 

Assume  c(j)  >-  0  ->  X(j)  starts  at  upper  bound 
inc  ( j)  =  -1 

if  (c(j)  -It.  zero)  inc(j)  =  1 
do  4505  i  -  1,  m 

AXINC (i»  j)  -  A(i, j) *DBLE(inc( j) ) 

( INC ( j)  -  -1) 

4505 

continue 

4510 

continue 

c 

c 

Now  calculate  RA,  the  matrix  of  ratios 
do  4520  i  =  1,  m 

do  4515  j  =  2,  n 
RA(i,j)  =  zero 

if  (A(i,j)  .ne.  zero)  RA(i,j)  -  A(i, j-1) /A(i, j) 

4515  continue 

if  (xprint(7)  .eq.  1)  then 

write  (8,*)  'ARATIO:  RA-> ' ,  (RA(i,  j)  ,  j«l,n) 
write  (8,*)  'ARATIO:  AXINC- ',  (AXINC (i, j) , j=l,  n) 
endif 

4520  continue 

if  (xprint(7)  .eq.  1) 

write  (8,*)  'ARATIO:  Inc-' ,  (inc  ( j) ,  j**l,n) 

return 

end 

c - 

c 

subroutine  CUTPT2 
c 

c  Called  by  PRECUT;  calc,  (alpha*)  -  int .  pts.  of  rays  &  I-cut 
c  See  paper  by  Balas,  et  al.  (BBGS71)  on  using  DUAL-UHC [XO] . 

c - 

include  ' $DISK2 : [SALT2 . ILPlJplist . for ' 
include  '$DISK2: [SALTZ . ILPlJcomprt . for • 
include  '$DISK2: [SALTZ . ILP1) comlpl . for ’ 
include  '$OISK2: (SALTZ. ILPl]comxrun2. for' 
c 

Integer  indr(nn) 

Real*8  f (nn) ,  fe2 (nn) ,  sigma (nn) , L(nn) , LIrO, vec (nn) , theta  (nn) 
Logical  nplus(nn) 
dn2  -  DBLE(n) /2.D0 
do  4600  j  =  1,  n 

f(j)  -  X0(j)  -  DINT  <X0 ( j ) ) 
fe2  { j)  =  f(j)  -  0.5D0 
indr(j)=  0 
4600  continue 

c 

c  Loop  through  all  of  the  extreme  directions 

do  4690  k  -  1,  nx 
NL  -  0 

do  4610  j  =  1,  n 
L ( j)  «  bigm 
nplus(j)  =  .false. 
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if  (zero  .It.  <fe2(j)*ABAL(j,k)))  then 

nplus(j)  -  -true. 

NL  -  NL  +  1 

L  ( j)  -  fe2(j)/ABAL(j,k) 
endif 

4610  continue 

c 

LIrO  -  zero 

if  (NL  .eq.  0)  goto  4640 
call  RSORT1 (n,L,  indr) 
c 

do  4630  i  =  1,  NL 
do  4625  j  -  1/  n 

vec(j)  -  fe2 ( j)  -  ABAL( j,k) *L(indr (i) ) 

4625  continue 

vecnrin  «  V1NORM  (n,  vec) 
if  (vecnrin  .gt.  dn2)  goto  4640 
LIrO  -  L(indr (i) ) 

4630  continue 

c 

4640  do  4650  j  -  1,  n 

theta (j)  -  f(j)  -  ABAL(j,k)*LIr0 

4650  continue 

c 

nnl  =  0 

do  4660  j  “  1/  n 
3igma(j)  -  -1.D0 
if  ( (theta ( j ) .gt . 0 . 5D0) . or . 

( (theta ( j) .eq.0.5D0) .and. (ABAL ( j , V) . It . zero) )) then 

sigma  (j)  =  one 
nnl  =  nnl  +  1 
endif 

4660  continue 

c 

tempi  =  VDOT(n, f, sigma) 
temp2  =  VDOT ( n , ABAL ( 1 , k ) , sigma ) 

LAMDA ( k ) =  (templ-dble (nnl) ) /temp2 
do  4670  j  =  1,  n 

ALPHA (j,k)  =  X0(j)  -  LAMDA (k)*ABAL(j,k) 

4670  continue 

if  (xprint(4)  .eq.  1)  then 

write (4,  *)  - - >  Ray  k  *  ’'k 

write  (4, *)  ’ABAL  ' ,  (ABAL ( j, k) , j-1, n) 
write(4,*)  'nplus  ' , (nplus ( j) » j=lrn) 
write  (4, *)  'L  ' »  (L( j) , j=lrn) 
write  (4, *)  'indr  • , (indr ( j) , j-l,n) 
write(4,*)  'LIrO  ',LIR0 
write (4, *)  'sigma  (sigma (j) , j=l»n) 
write(4,*)  'theta  ',  (theta (j) , j-1* n) 
write (4, *)  'nnl  ',nnl 
write  (4, *)  'LAMDA  '^LAMDA(k) 
write (4/ *)  'ALPHA  ' , (ALPHA( j,k) , j=l,n) 
endif 

4690  continue 
return 
end 
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c - 

c 

subroutine  REDCTS 
c 

c  Called  by  RUN  prior  to  XRUN. 

c  Eliminates  redundant  cts.  not  intersecting  n-simplex  S. 

C  THIS  ROUTINE  MODIFIES  PROBLEM  DATA. 

c - 

include  *$DISK2: [SALTZ . ILPl]PLIST.FOR' 
include  '$DISK2: [SALTZ.ILP1]C0MPRT.F0R' 
include  ' SDISK2 : [SALTZ . ILP1] comlpl . for 1 
include  'SDISK2: [SALTZ. ILP1] comrunl . for' 
include  *$DISK2: [SALTZ . ILP1] comxrun2 . for • 
include  *$DISK2: [SALTZ . ILP1] comxrun3 . for ’ 
include  '$DISK2: [SALTZ. ILP1] comxrun4 . for' 
logical*2  active (mm) 

if  (hprint (17) .eq.l)  write  (8,*)  'REDCTS:  eliminating  cts' 
c 

c  nactiv  *  counts  the  number  of  active  (non- redundant)  cts. 

nactiv  =  0 
c 

c  Check  each  constraint  for  redundancy 

do  6970  i  -  1,  m-2 
active (i)  -  .false, 
if  (BFCXO(i))  then 

c  Constraint  binding  at  XO:  always  active 

active(i)  «  .true, 
nactiv  -  nactiv  +  1 

if  (hprint (17) .eq.l)  write  (8,*)  '  ct.',i,'  active' 
goto  6970 
endif 
c 

c  Nonbinding  ct:  Does  it  cause  a  P(k)  to  be  infeasible? 

do  6960  k  »  1,  nx 
aipk  *  zero 
do  6957  j  -  1,  n 

aipk  «  aipk  +  A(i,  j) *P ( j, k) 

6957  continue 

c 

c  Test  feasibility  of  P(k)  wrt  this  constraint 

if  (aipk  .gt.  B(i))  then 

c  P  (k)  violates  (i) ,  so  (i)  is  active 

active (i)  =  .true, 
nactiv  -  nactiv  +  1 

if  (hprint (17) .eq. 1)  write  (8,*)  '  ct.',i, '  active' 
goto  6970 
endif 
c 

6960  continue 

c 

6970  continue 
c 

c  Make  sure  last  2  contraints  (I-Cut  t  OBJ)  are  active 

active (m-1)  «  .true, 
active (m)  -  .true, 

nactiv  -  nactiv  +  2 

c  If  all  original  cts.  binding  at  X0,  then  none  redundant 

if  (nactiv-2  .eq.  m)  goto  6990 
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c 

c  Delete/overwrite  redundant  eta.  from  XRUN  global  vara, 

ia  **  0 

do  6980  i  »  1,  ■ 

if  (.not.  active (i) )  goto  6980 
ia  =  ia  +  1 
B (ia)  -  B (i) 

INDCT(ia)  =  INDCT(i) 

BFCXO(ia)  -  BFCXO(i) 

AIMIN(ia)  -  AIMIN(i) 

AIMAX(ia)  =  AIMAX(i) 

PRICES (ia)=  PRICES (i) 


do  6975  j  =  1, 

n 

A(ia, j) 

-  A(i, j) 

c 

AXINC (ia, j ) 

-  AXINC  (i,  j) 

c 

RA(ia, j) 

-  RA(i,j) 

c 

CMPMIN (ia, j 

)  =  CMPMIN <i,  j) 

6975 

continue 

do  6977  k  -  1, 

nx 

CTXPT (k, ia) 

-  CTXPT (k,i) 

6977 

continue 

c 

CMPMIN (ia,N+l) 

=  CMPMIN (i,N+l 

6980 

continue 

c . 

6990 

M  =  nactiv 

c . 

return 

c - 

end 

C 

subroutine  REORDR 
c 

c  Called  once  per  search  hp  by  XCP  to  reorder  the  variables, 
c  Variables  fixed  at  a  value  (by  SHR(hp)  bounds)  go  fir3t. 
c  THIS  ROUTINE  MODIFIES  PROBLEM  DATA.  (See  Section  6.5) 

- - 

include  '$DISK2: [SALTZ . ILP1] plist . f or ' 
include  'SDISK2: [SALTZ . ILP1] comprt . for ' 
include  ' $DISK2 : [SALTZ . ILP1 ] comlpl . f or ' 
include  ' $DISK2 : [SALTZ . ILP1] comrunl . for ' 
include  '$DISK2: [SALTZ . ILP1] comxrun2 . for ' 
include  '$DISK2: [SALTZ . ILP1] comxrun3 . for ' 
include  ' $DISK2 : [SALTZ . ILP1] comxrun4 . for ' 

C 

logical*2  fxdvar(nn) 

real*8  tC(nn) ,tX0 (nn) ,tP (nn,nn) ,tDIRS ;nn,nn) 
real*8  t ABAL (nn, nn) , tALPHA (nn, nn) 

real*8  tA(mm,nn) , tLL (mm+1, nn) ,tUU(mm+l,nn) ,tCTXPT (nn,mm) 
integer  tSUB (nn) , tLONS (nn) , tUPNS (nn) , tLOBD (nn) , tUPBD (nn) 
integer  tXSTAR(nn) ,tINC (nn) 
c  tVAR  -  temporary  copy  of  variable  VAR 

if  (hprint (18) .eq. 1)  write  (8,*)  'REORDR:  reordering  vars 
c 

c  varord:  all  fixed  vars.  first,  followed  by  remaining  vars 

nfixed  =  0 
c 
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c  Check  all  the  variables  for  being  fixed 

do  6910  j  =  1,  n 

fxdvar(j)  -  .false, 
if  (lons(j)  .eq.  upns(j))  then 
fxdvar(j)  =  .true, 
nfixed  =  nfixed  +  1 
varord (nf ixed)  —  j 

if  (hprint (18) .eq. 1) write (8, *) '  Var ' , j, ' f ixed»* , Ions ( j ) 
endif 

6910  continue 

c 

c  If  nfixed  =  0,  then  variables  need  not  be  reordered 

if  (nfixed  .eq.  0)  goto  6945 
c 

c  Otherwise,  put  free  variables  after  fixed  variables 

jfree  =  0 
do  6915  j  =  1,  n 

if  (.not.  fxdvar(j))  then 
jfree  =  jfree  +  1 
varord(nf ixed+jf ree)  =  j 
endif 

6915  continue 

if  (hprint (18) .eq. 1)  write(8,*)'  varord  ' , (varord (L) , L«l, n) 
c 

c  Reorder  global  variables  into  temporary  global  vars. 

do  6930  j  =  1,  n 

jv  =  varord (j) 
tC ( j )  =  C ( jv) 

tX0(j)  =  X0(jv) 
tSUB(j)  =  SUB ( jv) 
tLONS(i)  =  LONS ( jv) 
tLOBD(j)  =  LOBD(jv) 
tUPNS(j)  =  UPNS(jv) 
tUPBD(j)  -  UPBD(jv) 
tXSTAR ( j ) =  IXSTAR(jv) 
tINC(j)  -  INC ( jv) 
c 

do  6920  k  =  1,  nx 

tP  ( j,k)  -  P ( jv,k) 

tALPHA(j,k)  -  ALPHA (jv,k) 
tDIRS ( j,k)  =  DIRS ( jv, k) 
tABAL ( j, k)  -  ABAL ( jv, k) 

6920  continue 

c 

do  6925  i  -  1,  m 

tA(i, j)  =  A (i, jv) 
tLL(i, j) =  LL(i,jv) 
tUU(i, j)=  UU(i, jv) 
tCTXPT ( j, i) =CTXPT ( jv, i) 

6925  continue 

tLL(m+l, j)«LL(m+l, jv) 
tUU(m+l, j)=UU(m+l, jv) 

6930  continue 
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Reorder  global  variables  from  temporary  global  vara 
do  6940  j  -  1,  n 


6932 

c 


6936 


c 

6940 

c 

c 


c 

6945 


C<j) 

X0<j) 

SUB ( j) 
LONS(j)  - 
LOBD(j)  - 
UPNS(j)  - 
UPBD(j)  - 
IXSTAR< j)- 
INC<j) 


tC(j) 
tXO(j) 
tSUB ( j ) 
tLONS ( j) 
tLOBD(j) 
tUPNS ( j) 
tUPBD  <  j ) 
tXSTARf j) 
tINC(j) 


do  6932  k  «  1,  nx 

P(j,k)  -  tP(j,k) 

ALPHA( j,k)  -  t ALPHA ( j,  k) 

DIRS ( j,  k)  -  tDIRS ( j, k) 

ABAL ( j ,  k )  -  tABAL<j,k) 

continue 


do  6936  i  «  1,  m 
A(i,j)  =  tA (i, j ) 

LL(i, j)«  tLL (i, j) 
UU(i,j)=  tUU (i, j) 

CTXPT ( j, i) -tCTXPT  <  j , i) 
continue 

LL (m+1, j )  -  tLL(m+l, j) 

UU (m+1, j)  -  tUU <m+l, j) 

continue 


Recalculate  other  global  variables  (CMPMIN,  INC,  RA 
call  MINCMP 
call  ARATIO 

continue 

return 

end 


i 


AXINC} 


61 


c - 

c 

subroutine  XRUN 
c 

c  Called  by  RUN  to  organize  Exact  Ceiling  Point  Algorithm  (XCPA) , 
c  given  initial  solution  &  value  {X*,Z*).  (See  Section  5.5,  Step  3) 
c 

c - 

include  '$DISK2: [SALTZ . ILPlJplist . for ' 
include  '$DISK2: [SALTZ . ILP1] comprt .for’ 
include  ’$DISK2: [SALTZ . ILP1] comlpl . for ’ 
include  *$DISK2: [SALTZ . ILP1] comrunl . for ’ 
include  '$DISK2: [SALTZ . ILP1] comxrun3 . for • 
logical*2  hplist (mm) 

if  (xprint (12) .eq. 1)  write(9,*)  'XRUN:  Z*  -  ',ZSTAR 
c 

icase  -  1 

c  Initialize  hplist  «  list  of  constraint  hp's  already  searched 

do  6001  i  “  1,  m 

hplist(i)  =  .false. 

6001  continue 

c 

c  Loop  thru  constraints,  enumerating  1-Ceiling  Points 

do  6080  i  »  1,  m 
c 

c  Select  a  search  hyperplane  (ihp) 

ihp  -  IXPICK (hplist) 

if  (xprint (12) .eq. 1)  write (9,*)  'XRUN:search  ct  “'/ihp 
hplist (ihp)  =  .true, 
if  (ihp  .eq.  0)  goto  6090 
c 

c  Enumerate  feasible  l-CP(ihp)'s  better  than  incumbent 

call  XCP(ihp) 
c 

if  (xprint (12) .eq. 1) write (9, *) 'xcalls ' , (xcalls (L) , L»l,n) 
if  (zstar  .ge.  zup)  goto  6090 
if  (icase  .eq.  2)  goto  6090 
c 

c  Reoptimize  (LPr) '  =  tightened  LP-relaxation  of  ILP 

call  XREOPT(ihp) 
c 

6080  continue 

c 

6090  continue 

return 
end 

c - 

c 

subroutine  XCP(hp) 
c 

c  Called  by  RUN  &  RUNCUT  to  find  ceiling  points  wrt  constraint  (hp) . 
c  (See  Section  5.4) 

c - 

include  * $DISK2 : [SALTZ. ILPlJplist. for’ 
include  'SDISK2: [SALTZ. lLPl]comprt. for' 
include  '$DISK2: [SALTZ . ILP 1] comlpl . for ' 
include  '$DISK2: [SALTZ . ILP1) comrunl . for ' 
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include  '$DISK2: [SALTZ . ILP1] comxrun2 . f or ' 
include  '$DISK2: [SALTZ . ILP1] comxrun3 . for ’ 
include  '$DISK2: [SALTZ . ILP1) comxrun4 .for’ 
integer  hp 
if  (xprint (13) .eq. 1) 

-  write (9, *)  'XCP:  Z*  -',ZSTAR, '  hp  =  ',hp 
c 

B(m)  »  -(zstar  +  all) 
do  6100  i  =  1,  m 
GAP ( i ,  1 )  =  B(i) 

6100  continue 

if  (xprint (13) .eq. 1) write  (9,*)  'gap (, 1) - ' , (gap (i, 1) , i«l,m) 
c 

iempty  =  ISHR(hp) 

if  (1  .eq.  iempty)  goto  6105 

if  (xprint (8)  .eq.OJwrite  (*,*)  'XCP:  NO  INT.  SOLNS.  in  SHR* 
if  (xprint (13) .eq. 1) write  (9,*)  'XCP:  NO  INT.  SOLNS.  in  SHR’ 
goto  6190 
c 

c  Put  the  fixed  vars .  first,  if  not  searching  I-Cut 

6105  if  (hp  .ne.  m-1)  call  REORDR 

c 

c  Initialize  key  arrays 

do  6110  j  -  1,  n 

upc(j)  =  upbd(j) 
loc(j)  =  lobd(j) 

c  upc ( j) ,loc( j)  -  upper  6  lower  conditional  bounds  on  X(j) 

first ( j)  =  0 
final(j)  *  0 

c  f ir3t  (j) , final (j)  -  interval  over  which  X(j)  will  range 

ISN(j)  -  0 

c  ISN  -  current  integer  solution  being  spelled  out 

newz(j)  =  0 

c  newz(j)  =  1  =>  New  Z*  found  with  previous  partial  solution 

xcalls(j)*  zero 

c  xcalls(j)=  counter  for  number  of  calls  to  XCB(j) 

6110  continue 

c 

lj  -  1 

if  (ICHKBD  <1 j)  -eq.  0)  goto  6190 
call  LOOPBD (1 j) 
c 

c .  Main  Enumeration  Loop  .  (See  Section  5.3)  . 

c 

6120  ISN(lj)  =  ISN(lj)  +  INC (1 j) 

c 

c  Take  a  forward  step  &  calc,  conditional  variable  bounds 

lj  -  lj  +  1 

call  XCB(lj) 

if  (ICHKBD  <1 j)  .eq.  0)  goto  6140 
c 

if  (lj  .It.  n)  then 

c  Set  the  loop  bounds  and  take  a  forward  step 

call  LOOPBD (lj) 
goto  6120 
endif 
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c 

c  Feasible  Completion 

6130  continue 

iempty  =  ILASTV(hp, 1 j) 
if  (iempty  .eq.  0)  goto  6190 
c 

c  Backtrack  to  lower  level 

6140  lj  -  lj  -  1 

if  (IBOTLP (1 j)  .eq.  1)  goto  6120 
if  (lj  .It.  2)  goto  6190 
goto  6140 
c 


c .  End  of  Main  Loop 

6190  continue 


c 

return 

end 

c - 

c 

subroutine  XCB(j) 
c 

c  Called  by  XCP  to  compute  bounds  for  X ( j )  conditioned  on  the  values 
c  of  ( X ( 1) ,  X(2),...,  X(j-l)).  Also  called  by  IBOTLP  after  a  new 
c  incumbent  has  been  found.  X(j)  corresponds  to  ISN(j) 
c  (See  Section  5.4.1) 

c - 

include  ' SDISK2 : [SALTZ . ILPlJplist . for ' 
include  ' SDISK2 : (SALTZ . ILP1) comprt . for ■ 
include  ' SDISK2 : (SALTZ . ILP1] comlpl . for ' 
include  ' SDISK2 : (SALTZ . ILP1] comxrun3 . for ' 
include  ' SDISK2 : [SALTZ . ILP1] comxrun4 . for ’ 
integer  j 
c 

if  (xprint(14)  .eq.  1)  write  (9,*)  'XCB:  j  =  ',j 
xcall3(j)  =  xcalls(j)  +  one 
c  if  (j  .eq.  1)  goto  6220 

if  (ISN(j-l)  .eq.  first(j-l))  goto  6205 
if  (newz(j)  .eq.  1)  goto  6205 
c 

c  Use  additive  form  after  first  value  of  ISN(j-l) 

do  6203  i  *  1,  m 

gap  (i, j)  =  gap ( i , j )  -  axinc(i,j-l) 
if  (a  (i, j) )  6201,  6203,  6202 

6201  LL  ( i ,  j )  ■=  LL  ( i ,  j )  +  RA(i,j) 
goto  6203 

6202  UU  (i,  j)  -  UU  (i,  j)  +  RA  (i,  j) 

6203  continue 
goto  6230 

c 

c  On  first  value  of  ISN(j-l),  use  multiplicative  form 

6205  do  6210  i  =  1,  m 

gap(i, j)  =  gap (i, j-1)  -  A(i, j-1) *ISN( j-1) 

6210  continue 

c 
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c  Compute  conditional  bounds  wrt  each  constraint  i 

6220  do  6225  i  =  1,  m 

if  (A(i, j) )  6221,  6225,  6222 

6221  LL(i,j)  “  (gap (i, j )  -  CMPMIN(i, j+1) ) /A(i, j) 
goto  6225 

6222  UU(i, j)  -  (gap  <i, j)  -  CMPMIN(i, j+1) > /A(i, j) 

6225  continue 

if  (xprint(14)  .eq.  1)  then 

write (9, *) 'xcb:LLj ' , (LL (i, j ) , i=l,m+l) 
write  (9,  *)  •xc^UUj',  (UU  (i,  j) ,  i=l,m+l) 
endif 
c 

c  Pick  the  tightest  condl .  bound  (row  M+l  < — >  uncdl.  bds) 


6230 

maxarg  =  m+l 
minarg  =  m+l 
do  6250  i  =  1, 

m 

if  (A(i, j) ) 

6235, 

6250,  6245 

6235 

if  (LL (i, j) 
goto  6250 

.gt. 

LL (maxarg, j) ) 

maxarg  = 

6245 

if  (UU (i,  j) 

.It. 

UU (minarg,  j) ) 

minarg  = 

6250 

continue 

upc(j)  =  IRNDWN (UU (minarg, j) ) 
loc(j)  =  IRNDUP (LL (maxarg, j) ) 
return 
end 

c - 

c 

Integer  Function  ICHKBD(j) 
c 

c  Called  by  XCP  to  check  if  ISN(j)  is  at  its  final  bound, 
c  Returns  1  if  loc(j)  <=  upc(j)  =>  take  forward  step 
c  0  "  "  >  "  =>  take  backward  step 

c - 

include  '$DISK2: [SALTZ . ILP1] plist . f or ' 
include  ' SDISK2 : [SALTZ . ILP1] comprt . for ' 
include  '$DISK2: [SALTZ . ILP1] comxrun3 . for ' 
include  '$DISK2: [SALTZ . ILP1) comxrun4 . for ' 
integer  j 
c 

c  if  (xprint(14)  .eq.  1) 

c  -  write  (9,*)  ' ICHKBD  j  = ' , j ,  ’  loc ' , loc ( j) , '  upc',upc(j) 

c 

ICHKBD  =  1 

if  (loc(j)  .gt.  upc(j))  then 
c 

c  Can  backtrack  to  next  value  of  ISN(j-l) 

ICHKBD  -  0 

if  (RA (minarg, j)  .gt.  RA (maxarg, j ) )  goto  6390 
c 

c  Can  double  backtrack  to  next  value  of  ISN(j-2) 

lj  =  j  -  1 
endif 
c 

6390  continue 
return 
end 
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c - 

c 

subroutine  LOOPBD(j) 
c 

c  Called  by  XCP  to  set  first  and  final  values  for  level  j. 

c - 

include  'SDISK2: [SALTZ . ILP1] plist . f or ' 
include  ' $DISK2 : [SALTZ . ILP1 ] comxrun4 . for ' 
integer  j 
c 

c  If  INC(j)  =  -1,  start  enumerating  from  upper  bound 

c  INC ( j )  is  set  in  ARATIO 

if  ( INC ( j )  .eq.  -1)  then 
first ( j)  =  upc ( j) 
final(j)  =  loc(j) 
c 

c  Otherwise,  start  enumerating  from  lower  bound 

else 

first ( j )  =  loc(j) 
final ( j)  =  upc(j) 
endif 

ISN ( j )  =  first (j)  -  INC ( j ) 

return 

end 

c - 

c 

Integer  Function  ILASTV(hp,j) 
c 

c  Called  by  XCP  when  a  feasible  completion  has  been  reached, 
c  Returns  0  if  SHR(ns)  or  SHR(hp)  is  empty. 

c - - - 

include  '$DISK2: [ SALTZ . ILP1 J plist . for ’ 
include  ' SDISK2 : [SALTZ . ILP1 ] comprt . for ' 
include  '$DISK2: [SALTZ . ILP1] comlpl . for ' 
include  '$DISK2: [SALTZ . ILP1 ] comxrun4 . for ■ 
integer  hp, j 

if  (xprint  (16)  .eq. 1)  write (9,*)  'LASTV:hp=' ,hp, '  j=',j 
c 

c  Assign  last  component  based  on  sign  of  C(n) 

ISN(n)  =  upc(n) 

if  <C(n)  .It.  zero)  ISN(n)  =  loc(n) 
do  6500  i  =  l,m 

gap(i,n+l)  =  gap(i,n)  -  A (i, n) *DBLE (ISN (n) ) 

6500  continue 

c 

c  Calculate  objective  function  value  of  new  incumbent 

z  =  zero 

do  6505  ij  =  1,  n 

z  =  z  +  C (i j) *DBLE (ISN(i j) ) 

6505  continue 

if  (xprint (16)  .eq.  1)  then 

write  <  9, *)  'ILASTV:  z=',z, '  IS= ' ,  (ISN (i j) , i j=l, n) 
write(9,*)  'ILASTV:  gap= ' , (gap (i, n+1 ) , i=l ,m) 
endif 
c 

c  Call  new  incumbent  routines 

ILASTV  =  INCMOD (hp, z) 
return 
end 
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- - 

c 

Integer  Function  IBOTLP(j) 
c 

c  Called  by  XCP.  Return  0  if  ISN(j)  has  reached  its  loop  limit 

c - - - 

include  *$DISK2: [SALTZ . ILPl]plist . for ’ 
include  '$DISK2: [SALTZ . ILP1] comxrun4 . for ' 
integer  j 
c 

IBOTLP  =  1 

if  (ISN(j)  .eq.  FINAL (j) )  then 
IBOTLP  =  0 
goto  6605 
else 

if  (newz(j)  .eq.  0)  goto  6605 
if  (j  .eq.  1)  goto  6605 

c  O/w,  New  Z*  just  found  and  j  .ne.  1,  so  get  new  conditional 

c  bounds  (See  Section  5.4.3). 

call  XCB(j) 
newz(j)  =  0 
endif 
c 

6605  continue 
return 
end 

- - 

c 

Integer  Function  INCMOD (hp, val) 
c 

c  Called  by  ILASTV  whenever  a  new  incumbent  is  found  (fairly  rare) . 
c  Returns  0  if  new  (possibly  smaller)  SHR  is  empty;  1,  otherwise. 

c - 

include  '5DISK2: [SALTZ . ILPl]plist . for ' 
include  'SDISK2: [SALTZ . ILP1] comprt . for ' 
include  ' SDISK2 : [SALTZ . ILP1J comlpl . for  * 
include  ' $DISK2 : [ SALTZ . ILP1] comrunl . for ' 
include  ' SDISK2 : [SALTZ . ILP1] comxrun4 .for' 
integer  hp,iret4 
real*8  val 

if  (xprint (16) .eq.ljwrite (9, *)  ' INCMOD :hp hp, '  val  ',val 

c 

iret4  =  1 

delz  =  val  -  ZSTAR 
if  (delz  .It.  tol)  goto  6790 
c 

c  We  have  an  improved  feasible  solution 

ZSTAR  =  val 
do  6700  j  =  1,  n 

ixstar(j)=  ISN(j) 
newz(j)  =  1 

gap (m, j)  =  gap (m, j)  -  delz 
6700  continue 

gap(m,n+l)  =  gap(m,n+l)  -  delz 

if  (xprint(8)  .eq.0)  write(*,*)  '  New  Z*  **',  ZSTAR 
if  (xprint (16) .eq.l)  then 

write (9,*)  '  Ne\.  Z*  -  ',zstar 
write(9,*)  '  New  X*  =  ' , (ixstar ( j) , j=l, n) 
endif 
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c 

B(m)  =  -(zstar  +  all) 
call  OBJCUT (zstar  +  all) 
iret4  =  ISHR(hp) 

C 

6790  INCMOD  *  iret4 
return 
end 

c - 

c 

integer  function  IXPICK(list) 
c 

c  Called  by  XRUN  to  select  search  constraint  hyperplane, 
c  Do  not  cut  any  deeper  than  CX  =  Z*  hyperplane.  (See  Section  5.7) 

c - 

c$include :plist . for 
c$include : comprt . for 
c$include : comlpl . for 
c$include : comrunl . for 

include  1 SDISK2 : (SALTZ . 1LP1 ] plist . for 1 
include  ' $DISK2 : [ SALTZ . ILP1 ] comprt . for ' 
include  1 SDISK2 : (SALTZ . ILP1] comlpl . for • 
include  '$DISK2: [SALTZ . ILP1] comrunl . for ' 
logical*2  list (mm) 
c 

c  li3t(i)  =  .true.  =>  hp  (i)  has  already  been  a  search  ct . 

if  (xprint(17)  .eq.l)  write  (9,*)  "IXPICK:  zup=  " , zup 
c 

IXPICK  =  0 
zbarmn  =  bigm 
c 

c  Find  ct .  (IXPICK)  yielding  largest  decrease  in  zbar' 

do  6810  i  =  1,  tn 

if  (.not.  BFCXO(i))  goto  6810 
if  (list(i))  goto  6810 

Opt.  price  on  ct.  i  comes  from  final  tableau. 

(Works  fine  if  ALL  cts.  are  .GE.  or  ALL  are  .LE.) 
delb  =  one  +  AIMAX(i) 
delzO  =  DABS (PRICES (i))*delb 
zbar  =  zO  -  delzO 
c 

if  (zbar  .It.  zbarmn)  then 

c  Have  found  a  new  lowest  upper  bound  zbar:' 

IXPICK  =  i 
zbarmn  =  zbar 
endif 

6810  continue 
c 

if  (xprint(17)  .eq.  1) 

write  (9,*)  'IXPICK:  prices*' , (prices (i) , i-l,m) 
if  (IXPICK  .eq.  0)  GOTO  6840 
C 

c  Assume  old  case:  decrease  zbar  as  much  as  possible 

zup  •  DINT (zbarmn) 

c  Check  for  new  case:  only  decrease  upper  bound  to  Z* 

if  (zup  .gt.  zstar)  goto  6820 
zup  *  zstar 

AIMAX (IXPICK)  =  (Z0  -  zstar) /DABS (prices (IXPICK) ) 
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c 

c  Replace  I-cut  with  ceiling  point  constraint  < — >  IXPICK 

6820  B(m-l)  =  B(IXPICK)  -  AIMAX(IXPICK) 

do  6830  j  =  1,  n 

A(m-l,j)  =  -A (IXPICK, j) 

6830  continue 

AIMAX(m-l)  =  AIMAX( IXPICK) 
c 

6840  if  (xprint(17)  .eq.  1) 

write  (9,*)  'IXPICK  ct:  IXPICK, '  =>  zup=',zup, '  Z*=',Z3tar 

return 

end 

c - 

c 

subroutine  XREOPT(hp) 
c 

c  Called  by  XRUN  to  reoptimize  (LPr) '  after  searching  SHR(hp) . 
c  (See  Section  5.7) 

- - 

include  ’$DISK2: [SALTZ . ILP1] plist . for ' 
include  '$DISK2: [SALTZ . ILP1 ) comprt . for  1 
include  '$DISK2: [SALTZ . ILP1] comlpl . for ' 
include  '$DISK2: [SALTZ . ILP1] comrunl . for ' 
integer  hp 

if  (xprint (8) .eq. 0)  write(*,*)‘  XREOPT:  hp=',hp 
c 

c  Translate  constraint  hyperplane  (hp)  by  t 

M  =  M  -  2 

t  =  AIMAX(hp)  +  one 
B  (hp)  -  B (hp)  -  t 
c 

c  Find  solution  to  (LPr)'  as  defined  in  Section  5.7 

c  XBAR'  =  XBAR  -  t*BINVERSE ( , hp) 

c  ZBAR '  =  ZBAR  -  t*PI(hp) 

C 

do  6850  i  =  1,  m 

XBAR<ibasis <i) )  =  XBAR (ibasis (i) ) -t *ABAR (i, INITBV (hp) ) 

6850  continue 

do  6855  j  =  1,  n 
X0(j)  =  XBAR(j) 

6855  continue 

Z0  =  Z0  -  t*PRICES (hp) 
c 

c  Do  some  PRECUT-type  operations  on  the  data 

call  OBJCUT (ZSTAR+one) 

C  Reappend  ceiling  point  condition  &  objective  fn.  constraints 

do  6860  j  =  1,  n 
A(m+1,  j)  =  zero 
A(m+2,  j)  =  — C  ( j ) 

6860  continue 

B(m+1)  =  zero 

B(m+2)  =  zero 

M  =  M  +  2 
ZUP  =  DINT (Z0) 

if  (xprint (8) .eq.0)  write(*,*)'  XREOPT:  Z0=',Z0 
C 

return 

end 
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c - - 

c 

function  ISHR(hp) 
c 

c  Called  by  XCP  to  calc,  unconditional  bounds  -  Search  HyperRectangle 
c  Returns  0  if  SHR  is  clearly  empty,  1  if  SHR  is  not  necc.  empty. 

c - 

include  '$DISK2: [SALTZ . ILPl]plist . for ' 
include  '$DISK2: [SALTZ . ILP1] comprt . for ’ 
include  ' SDISK2 : [SALTZ . ILPl]comlpl . for ' 
include  'SDISK2: [SALTZ . ILP1] comrunl . for ' 
integer  hp,iretl 

if  <xprint(9)  .eq.  1)  write<8,*)  'ISHR:  hp  -',hp 
c 

if  (hp  .ne.  m-1)  goto  5005 
c 

c  CUTSHR: special  SHR  routine  for  Search  I-Cut 

call  CUTSHR 
c 

iretl  =  ISREND(hp) 
goto  5010 
c 

c  Search  ct .  is  a  functional  ct .  or  objective  fn.  hp 

5005  iretl  =  ISHRSC(hp) 

5010  ISHR  =  iretl 

return 
end 

c - 

c 

integer  function  ISHRSC(hp) 
c 

c  Called  by  ISHR  to  calc,  uncondl.bds.  for  functional  search  ct.(hp). 
c  Returns  0  if  SHR  is  clearly  empty,  1  if  SHR  is  not  necc.  empty, 
c  (See  Section  5.3) 

c - 

include  '$DISK2: [SALTZ. ILPljplist .for' 
include  '$DISK2: [SALTZ . ILP1] comprt . for ' 
include  '$DISK2: [SALTZ . ILP1] comlpl . for ' 
include  ' SDISK2 : [SALTZ . ILP1] comrunl . for ' 
include  ' SDISK2 : [SALTZ . ILP1] comxrun2 . for ' 
include  '$DISK2: [SALTZ . ILP1] comxrun3 . for ' 
include  'SDISK2: [SALTZ . ILPl]comxrun4 .for ' 
integer  hp,iret2 

real *8  pts (nn, 3*nn) , a row (nn) , rmin, max 
if  (xprint(9)  .eq.  1) 

write  (8,*)  'ISHRSC:  hp,icase  =',hp,icase 
c 

c  First,  get  the  bounds  for  SHR (n-simplex) 

iret2  =  ISRSIM(m) 
ti  =  aimax(hp) 
c 

c  Store  coefficients  of  the  search  ct .  hyperplane 

do  5050  j  =  1,  n 

arow  ( j )  =■  A  (hp,  j) 

5050  continue 

c 
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c 


5055 

5060 

c 

c 


\ 


5062 


5063 

c 

c 


5070 

c 

c 


5080 

c 

c 


5082 

c 


nrow  =  0 

Get  points  where  kth  extreme  dir.  meets  ceiling  point  ct.(i') 
do  5060  k  =  1,  nx 

if  (CTXPT (k, hp) )  goto  5060 
nrow  =  nrow  +  1 
t  =  bigm 

tdot  »  VDOT (n,arow, dirs (1,  k) ) 
if  (tdot  .gt.  .001)  t  =  ti/tdot 
do  5055  j  =  1,  n 

pts(j,nrow)  =  X0(j)  -  t*dirs(j,k) 
continue 
continue 

Add  the  alpha (k)'s  and  P(k)'s  to  the  list  of  points 
do  5063  k  =  1,  n 

if  (.not.  CTXPT (k,hp))  goto  5063 
do  5062  j  -  1,  n 

pts ( j,nrow+l)  =  alpha  (j,k) 
pts ( j,nrow+2)  =  P(j,k) 
continue 
nrow  =  nrow  +  2 
continue 

For  each  component  j,  find  the  min  &  max  over  all  pts(,k) 
do  5080  j  =  1/  n 
rmin  =  pts(j,l) 
rmax  =  pts ( j, 1) 
do  5070  k  =  2,  nrow 

if  <pts(j,k)  .It.  rmin)  rmin  =  pts(j,k) 
if  (pts(j,k)  .gt.  rmax)  rmax  «  pts(j,k) 
continue 

Round  these  reals  to  integer  bounds 
lobd(j)  =  IRNDUP (rmin) 
upbd(j)  =  IRNDWN (rmax) 

Take  intersection  of  n-simplex  and  uncdl .  bds . 
lobd(j)  =  MAX0 (lobd( j) ,  Ions  ( j) ) 
upbd(j)  -  MIN0 (upbd ( j) ,  upns(j)) 
continue 

Check  for  icase  =  2  =>  (shr(ns))  =  (shr(hp)) 
if  (icase  .eq.  2)  goto  5088 
do  5082  j  =  1,  n 

if  ( (upns ( j ). eq . upbd ( j) ). and. (Ions ( j )). eq. lobd ( j) )  goto  5082 
icase  =  3 

if  (xprint(8)  .eq.  0)  then 

write(8,*)  '  ***  Case  C:  lshr(hp)i  <  |shr(ns)|  ***' 
write  (,*,*)  '  ***  Case  C:  |shr(hp)|  <  |shr(ns)|  ***' 
endif 
goto  5088 
continue 
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icase  =  2 

c  In  case  2,  remove  ceiling  point  constraint  ( i 1 ) 

B(m-l)  =  zero 
do  5084  j  =  1,  n 
A(m-l,j)  =  zero 
5084  continue 

AIMAX(m-l)  =  zero 
c 

5088  iret2  =  ISREND(hp) 

do  5090  j  =  1,  n 

irange(j)  =  irange ( j ) *iret2 
5090  continue 

c 

ISHRSC  =  iret2 

return 

end 

c - 

c 

integer  function  ISREND(hp) 
c 

c  Called  by  ISHR  to  complete  calc,  uncondl.bds. 
c  Returns  1  =>  SHR  non-empty;  =  0  =>  SHR  empty 

c - 

include  '$DISK2: [SALT2 . ILPllplist . for • 
include  ' SDISK2 : [SALTZ . ILPl] comprt . for ' 
include  ' $DISK2 : [ SALTZ . ILPl ] comlpl . f or  * 
include  ' SDISK2 : [SALTZ . ILPl] comxrun2 . for ' 
include  ' $DISK2 : [SALTZ . ILPl] comxrun3 . for ' 
integer  hp 
c 

c  Find  range (j)  :=  [U ( j) -L ( j) ] +1;  size  product  of  Ranges 

size  =  one 
do  5100  j  =  1,  n 

c  Prevent  uncondl .  bds .  from  exceeding  [0,  SUB] 

lobd(j)  =  MAX0 (lobd( j) ,  0) 
upbd(j)  =  MIN0 (upbd( j) ,  SUB  ( j) ) 
irange (j)  =  1+ (upbd ( j) -lobd ( j) ) 
size  =  size*DBLE (irange ( j) ) 

5100  continue 

if  (xprint(8)  .eq.  0) 

write(*,*)  *  ISREND:  3ize  =',size,'  sizlim  =',sizlim 
c 

if  (xprint(10)  .eq.  1)  then 

write  (8,*)  '  - ' 

write  (8,*)  '  Search  hp  '  ,hp 
write  (8,*)  '  Pts  in  hp  ',size 
do  5103  j  =  1,  n 

write  (8,*)  j, lobd (j) , upbd (j) 

5103  continue 

write  (8,*)  '  - ' 

endif 

c 

c  Exit  if  the  search  hyperrectangle  is  empty 

ISREND  =  0 

if  (size  .le.  zero)  goto  5140 
c 
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c  If  it  ia  non-empty,  (re) initialize  key  variables 

ISREND  =  1 
if  (hp  .It.  m)  then 
c  Compute  CMPMINU,  j) 

call  MINCMP 
c 

c  Initialize  Conditional  variable  bounds  wrt  each  ct . 

do  5120  i  =  1,  m+1 
do  5115  j  =  1,  n 

LL(i, j)  =  DBLE  (LOBD ( j ) ) 

UU  (i, j)  =  DBLE (UPBD  ( j) ) 

5115  continue 

5120  continue 

if  (size  .It.  sizlim)  goto  5140 

STOP  '  ****  Search  Rectangle  size  exceeds  sizlim  ***♦' 

endif 

5140  continue 
return 
end 

c - 

c 

integer  function  ISRSIM(hp) 
c 

c  Called  by  ISHRSC  and  CUTSHR  to  calculate  unconditional  bounds, 
c  Returns  1  =>  SHR  non-empty;  0  =>  SHR  empty. 

c - 

include  '$DISK2: (SALTZ . ILP1] plist . f or ' 
include  'SDISK2: [SALTZ . ILP1J comprt . for ' 
include  '$DISK2: [SALTZ . ILP1] comlpl . for ' 
include  '$DISK2: [SALTZ . ILP1] comxrun2 . for • 
include  'SDISK2: [SALTZ . ILP1] comxrun3 . for ' 
integer  hp,iret3 

real*8  shr (2 *nn, nn) , shrmin, shrmax 

if  (xprint(lO)  .eq.  1)  write  (8,*)  'ISRSIM:  hp=',hp 
c 

c  Append  P(j,k)  to  shr,  the  matrix  of  points  used  to  form  SHR 

nrow  =  0 

do  5160  k  =  1,  nx 

c  if  hp  =  m,  all  ad j .  ext.  pts.  are  used  to  form  SHR(ns) 

if  (hp  .eq.  m)  goto  5148 
if  (.not.  ctxpt(k,hp))  goto  5160 
5148  nrow  =  nrow  +  1 

do  5150  j  =  1,  n 

shr(nrow,j)  =  P(j,k) 

5150  continue 

c 

c  Append  ALPHA (k)  to  shr  if  it's  closer  to  X0  than  P(k) 

if  (DABS (X0 ( 1 )  -  ALPHA ( 1, k) )  .It. 

DABS (X0 ( 1)  -  P(l,k)))  then 
nrow  =  nrow  +  1 
do  5155  j  =  1,  n 

shrfnrow, j)  =  ALPHA ( j , k) 

5155  continue 

endif 

5160  continue 
c 
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c  Now  get  bounds  from  SHR 

do  5170  j  -  1,  n 

shrmin  =  shr(l,j) 
shrmax  =  shr (1, j) 
do  5165  i  =  2,  nrow 

if  (shr(i,j)  .It.  shrmin)  shrmin  =  shr(i,j) 
if  (shr(i, j)  .gt.  shrmax)  shrmax  =  shr(i,j) 

5165  continue 

c 

lobd<j)  =  IRNDUP (shrmin) 
upbd(j)  =  IRNDWN (shrmax) 

5170  continue 

c 

iret3  =  ISREND(hp) 
do  5175  j  =  1,  n 

irange(j)  =  irange ( j ) *iret3 

5175  continue 

c 

ISRS1M  =  iret3 
if  (hp  .eq.  m)  then 
do  5180  j  =  1,  n 
lons(j)  =  lobd(j) 
upns ( j )  =  upbd ( j ) 

5180  continue 

endif 
c 

return 

end 

c - 

c 

subroutine  MINCMP 
c 

c  Called  by  ISREND  &  REORDR  to  calc.  CMPMIN  =  Matrix  of  min.  completions 
c  CMPMIN ( i, j ) =  SUM  from  k=j+l  to  n  of  Min{ A (i, k) *UPBD (k) , A (i, k) *LOBD (k) ) 
c  (=  ”w(i,j)"  in  Section  5.4) 

c - 

include  '$DISK2: [SALTZ . ILPl]plist . for ’ 
include  '$DISK2: [SALTZ . ILPl]comprt . for ' 
include  ' SDISK2 : (SALTZ . ILP1] comlpl . for ' 
include  '$DISK2: [SALTZ . ILP1] comxrun2 . for ' 
include  '$DISK2: [SALTZ . ILP1 ] comxrun3 . for ' 
c 

do  5220  i  =  1,  m 
c 

c  (Initialize  CMPMIN  for  the  multirun  case) 

CMPMIN (i,n+l)  =  zero 
c 

c  More  efficient  to  start  from  last  column  &  work  toward  first 

do  5215  j  =  n,  1,  -1 

if  (A(i, j) )  5205,  5205,  5210 
c  if  A(i, j)  <  or  =  0 : 

5205  CMPMIN (i, j)  =  CMPMIN (i, j+1)  +  A (i, j) *DBLE (upbd ( j) ) 

goto  5215 

c  if  A(i, j)  >  0 : 

5210  CMPMIN (i, j)  -  CMPMIN (i, j+1)  +  A(i, j) *DBLE(lobd( j) ) 

5215  continue 

if  (xprint (10) .eq. 1) write (8, *) '  CMPMIN' , (cmpmin (i, j) , j=l,n+l) 
5220  continue 
return 
end 
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SWITCHES. U AT 


0000000000000000001010000 
0000000100000000000000000 
10.D16  1 .DIO 


1  2 

34567890123456 

7  8  9  0 

1  2  3  4  5 

1st 

row  is  hprint : 

2nd 

row  is  xprint: 

i 

1. 

getabc  (data  echo) 

1. 

Runcut 

2. 

z2phase 

2. 

Precut 

I 

3. 

zsolve 

3. 

Ob j cut,  Redcts 

4. 

zsetrs,  zpivot4 

4. 

Cutpt,Cutpt2 

i. 

5. 

setup 

5. 

Cuthp 

6. 

balas 

6. 

Cutshr 

f 

7. 

hrun 

7. 

Aratio 

8. 

hsdir 

8. 

1->MINIMAL  SCREEN  PRINTING 

1 

9. 

findfs 

9. 

Ishr 

10. 

f is2hp 

10. 

Ishrsc,  Isrend,Mincir.p 

T 

11. 

if feas 

11. 

Bounds 

12. 

phase3 

12. 

Xrun 

13. 

stayf s 

13. 

Xcp 

14. 

twovar,  getfes 

14. 

Xcb 

15. 

hround,  hrndpt 

15. 

Ichkbd 

16. 

VarOrder:  l)lo-hi;  2)hi-lo 

16. 

Ilastv,  Incmod 

17. 

REDCTS :  list  active  cts. 

17. 

Ixpick 

18. 

REORDR:  fixed  vars.  first 

18. 

1  =>  STOP  AFTER  HRUN 

k 

19. 

HRUN:  time/print  Phases 

19. 

1  =>  STOP  AFTER  ICUT 

20. 

HRUN:  1  =*>  nhps  =  No.BFCXO 

20. 

1  =>  only  use  cutpt2 

21. 

HRUN:  1  =*>  avoid  NCP  check 

21. 

Skip  HRUN,  start  RUNCUT  w/ 

22. 

22. 

23. 

23. 

24. 

24. 

25. 

25. 

3rd 

row:  sizlim  (for  SHR  size), 

callim 

(for  XCALLS)—  NOT  USED 

Note:  for  RANDOMLY  GENERATED  PROBLEMS:  set  hprint(16)  -  2 
for  REALISTIC  "  :  -  0 


ILFDATA, EAT  (sample  file) 

FC10 -DAT 
END 
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Example:  ("FC-IO"  from  Trauth  and  Woolsey,  1969) 
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Abstract 

An  Exact  Ceiling  Point  Algorithm 
for  General  Integer  Linear  Programming 

Robert  M.  Saltzman  and  FYederick  S.  Hillier 
Stanford  University,  1988 

This  report  describes  an  exact  algorithm  for  the  pure,  general  integer  linear  program¬ 
ming  problem  ( ILP ).  Common  applications  of  this  model  occur  in  capital  budgeting 
(project  selection),  resource  allocation  and  fixed-charge  (plant  location)  problems.  The 
central  theme  of  our  algorithm  is  to  enumerate  a  subset  of  all  solutions  called  “feasible 
1-ceiling  points.”  A  feasible  1-ceiling  point  may  be  thought  of  as  an  integer  solution  lying 
on  or  near  the  boundary  of  the  feasible  region  for  the  LP-relaxation  associated  with  (ILP). 
Precise  definitions  of  1-ceiling  points  and  the  role  they  play  in  an  integer  linear  program 
are  presented  in  a  recent  report  by  the  authors.  One  key  theorem  therein  demonstrates 
that  all  optimal  solutions  for  an  (ILP)  whose  feasible  region  is  non-empty  and  bounded 
are  feasible  1-ceiling  points.  Consequently,  such  a  problem  may  be  solved  by  enumerating 
just  its  feasible  1-ceiling  points.  Our  approach  is  to  implicitly  enumerate  1-ceiling  points 
with  respect  to  one  constraint  at  a  time  while  simultaneously  considering  feasibility.  Com¬ 
putational  results  from  applying  this  incumbent-improving  Exact  Ceiling  Point  Algorithm 
to  48  test  problems  taken  from  the  literature  indicate  that  this  enumeration  scheme  may 
hold  potential  as  a  practical  approach  for  solving  problems  with  certain  types  of  structure. 
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